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PREFACE 


The principal aim of this monograph i6 to present a coherent and fairly 
comprehensive account of recent progress in the area of dynamics and stability 
of mechanical systems with follower forces. By "recent, 11 quite specifically, 
is meant the period after 1963, the year of publication of the English trans- 
lation of the first book (by V. V. Bolotin) devoted in its entirety to non- 
conservative problems of the theory of elastic stability, i.e., problems with 
follower forces. 

The last decade has witnessed a considerable expansion of interest in 
this problem area, but the progress has been reported piecemeal by a variety 
of investigators in different countries and scattered in numerous journals. 

Even though advances are being continually made, it still appears to be 
justified to attempt to present an account of recent developments and to 
place them into a relative perspective. In this attempt, the author's own 
work and that of his collaborators has received, quite naturally, particular 
emphasis. 

It is hoped that the monograph may prove useful as a source of information 
on the current state-of-the-art for the research worker and practicing engineer. 
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CHAPTER I 


INTRODUCTION 


1.1 Structural Stability; Column Buckling 

The structural engineer and the applied mechanician are usually becoming 
acquainted with the area of structural stability through Euler's problem of 
elastic column buckling. 

There are several different ways in which the problem of column buckling 
can be presented to the beginner, but one of the most instructive ones is 
through an eccentrically loaded cantilevered column (Fig. 1.1) as was done 
recently by Ziegler [1].* 

It is assumed that the column is homogeneous, obeys linear elastic Hooke's 
law (Young's modulus E) and is of constant cross-section. Let i be the length 
of the column, El the flexural rigidity, e the eccentricity of the (compres- 
sive) load P and f the deflection at the free end. The load P is assumed to 
be of the "dead" type, i.e., a weight which does not change in magnitude and 
direction as a possible consequence of the deformation of the column. If the 
slope of the deflected column axis w(x) is assumed to be small as compared to 
unity, the bending moment at section x is 

M = P(e + f - w) (1.1) 


On the other hand in elementary theory of beam bending the bending moment M 
is known to be proportional to the curvature, i.e.. 


M = El w" 


( 1 . 2 ) 


Elimination of M in the above two relations leads to the differential equation 
of the deflection curve 

EI w # = P(e + f - w) (1.3) 

or 

w" + |r w = ~ (e + f) (1.4) 


Numbers in brackets indicate references compiled in a listing beginning 
on page 143. 



The general solution of this equation can be written in the form, with the 

o 

abbreviation P/(EI) * k 

w = A cos k x + B sin n x + e + f (1*5) 

The unknown constants A, B and f are to be determined from the boundary con- 
ditions 


w(0) * w'(0) “ 0 and w(£) = f 
which leads to the solution 


w 


cos kj£ 


(1- COS K x) 


and the end deflection 

f = e ( — - — - - 1) 

\ cos H A / 

It is seen from this last expression that if 


cos k ii c 1 ; i.e. h = » 


tt 2 EI 
4 A 2 


( 1 . 6 ) 


(1.7) 


( 1 . 8 ) 


(1.9) 


the deflections at the free end f become infinite, regardless how small the 
eccentricity e. The column is said to buckle under the critical load 
(buckling load) 


P 


1 


tt 2 EI 


4 tr 


( 1 . 10 ) 


Due to the assumptions introduced, the above relations are obviously 
valid only for small deflections. If a centrally loaded column is considered 
and if the analysis is based on a more exact (nonlinear) differential equa- 
tion of the deflection curve, which allows for large slopes of this curve, 
the dependence of P on the end deflection f can be established, with the re- 
sult illustrated in Fig. 1.2. 

In the range of the load 0 < P < P^ there is only one equilibrium posi- 
tion possible, namely, that of the straight column (f =0). This equilibrium 
position is stable in the sense that small disturbances or imperfections of 
various sorts are not followed by large deviations from this position. In 
the range P > P^, by contrast, the perfectly straight column (f = 0) is still 
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in equilibrium, but this equilibrium is unstable since small disturbances will 
cause the column to move away from this position (f # 0) . The stable equilib- 
rium positions in this range, P > P^, are located on a symmetric curve which 

branches from the straight line f * 0 at the point P « P^. 

It is thus seen that buckling of a column is associated with the pheno- 
menon of bifurcation of equilibrium paths, a concept intimately associated 
with Euler's notion of stability and instability. 

This concept of Euler in analyzing stability served technology well, par- 
ticularly in the area of structural engineering and structural mechanics, as 
applied to buckling of beams, frames, plates, etc., and various combinations 
of structural elements, subjected to dead loads. 

It was found, however, that this concept cannot be applied indiscrimi- 
nately to the stability problem of any mechanical system. Specifically, 

systems which are not subjected to dead loads but rather to forces due to an 

interacting medium have often to be analyzed differently with regard to sta- 
bility. Examples of such mechanical systems include airfoils placed in an 
airstream, turbine blades interacting with water, flexible pipes conveying 
fluid, elastic systems subjected to impinging fluid jets, as well as certain 
types of electro -mechanical interactions. 

A common feature of such mechanical systems, or rather of forces acting 
on them, is that these forces are not derivable from a potential (by contrast 
to dead loads) and as a result a stability analysis based on Euler's concept 
of bifurcation of equilibrium may break down. 

For such problems a more fundamental approach to problems of stability 
has to be followed, one which is, for instance, based on the dynamic method 
of investigating small motions induced as a result of perturbations of the 
positions of equilibrium. One finds then that stability of a mechanical 
system subjected to forces which are not derivable from a potential may, gen- 
erally, be lost either through oscillations with increasing amplitude or 
through a nonoscillatory motion away from the equilibrium position. In the 
former case no bifurcation of equilibrium exists and therefore Euler's method 
breaks down. In the latter case the dynamic method leads essentially to the 
same results as Euler's static method of analyzing stability. 

A significant property of forces which are not derivable from a potential 
is their dependence on the instantaneous position or configuration of the sys- 
tem upon which they are acting. That is these forces follow in some a priori 
prescribed manner the motion of the system. For this reason they have been 
termed not only nonconservative forces, but also circulatory forces, configu- 
ration-dependent forces or simply follower forces. 
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1.2 Aim and Scope of the Monograph 


The present monograph centers on problems of stability of equilibrium of 
mechanical systems with follower forces. Follower forces, as acting on me- 
chanical systems, may be of aerodynamic, hydrodynamic, electromagnetic or 
thermal origin. Furthermore, they occur frequently in automatic control 
systems. 

The beginnings of analyses of stability of mechanical systems with fol- 
lower forces go back to the late nine teen- twenties and are associated with the 
name of Nikolai [2,3] in Russia. Comprehensive, fundamental studies were car- 
ried out by Ziegler [4-7] in the fifties in Switzerland. The book by Bolotin 
[8], devoted in its entirety to nonconservative problems of the theory of 
elastic stability, presents a well-rounded state of knowledge as of a decade 
ago. 


Several areas of stability problems of mechanical systems with follower 
forces, such as the highly developed area of aeroelasticity (cf. Garrick [9]) 
and stability of rotating shafts, will not be considered in the present mono- 
graph since these areas have already received considerable attention. 

The primary purpose of the present monograph is limited in the sense that 
attention is confined to the developments of the last decade, i.e., after the 
publication of [8], and narrowed down further by emphasizing the analytical and 
experimental investigations in which the author and his coworkers were involved 
during the period of the last seven years. A review of the work, including 
numerous references, through the year 1966, is contained in reference [10]. 

Concepts of stability in mathematical terms, as well as criteria of sta- 
bility are reviewed briefly in the still introductory Chapter II, together with 
means of analytical specification of follower forces. Chapter III is devoted 
to a discussion of nondissipative (i.e., purely elastic) systems with two 
degrees of freedom. An illustrative example is considered first and a general 
linear system next. A remarkable feature of systems with follower forces is 
that even small damping forces and certain other velocity-dependent forces may 
have a strong destabilizing effect. Such destabilizing effects, both in dis- 
crete and continuous systems, are treated in detail in Chapter IV. The special 
considerations which must be introduced in the analysis of continuous systems 
are discussed in Chapter V. Mechanical systems with follower forces may re- 
quire particular procedures in their stability analysis. Such methods, in- 
cluding energy considerations, are dealt with in Chapter VI. 

The analytical work on systems with follower forces is sometimes being 
criticized as being purely mathematical and as having no relevance to actual 
mechanical devices and structures. To counter this argument, several possibi- 
lities of physical realization of mechanical systems with follower forces are 
examined in Chapter VII. Qualitative observations on demons trational labora- 
tory models and quantitative experiments are reported in Chapter VIII. 
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CHAPTER II 


CONCEPTS OF STABILITY AND FOLLOWER FORCES 


The term “stability* 1 assigns a quality to a state of a system which 
signifies that possible disturbances of the system will not essentially change 
the state. This qualitative description is necessarily vague and precise 
mathematical meaning is to be assigned to the terms “state, 11 “disturbances" 
and “essential change* “ 

The required mathematical apparatus has been supplied by Liapunov [11]. 

Let us consider a discrete system with r degrees of freedom described by r 
generalized coordinates q^, and let us examine the special case of the state 
of equilibrium 

q ± - 0 (i - l,2,*..,r) (2.1) 

If the system is disturbed at a time t =* t^, at any instant t its state will 
be characterized by (generally nonvanishing) coordinates q, and by generalized 
dq i 

velocities q^ = and can be thought of as a point in a 2r-dimensional 

Euclidean space with coordinates z^ 

z k " z k (t) (k “ 1 ’ 2 ’-*- 2r > < 2 - 2 > 


The state of equilibrium (2.1), according to Liapunov, Is said to be stable if 
for any e > 0 we can find a b > 0, depending on e only (and possibly on t Q ) 
such that 

2r 

Y z 2 < 6 at t = t Q (2.3) 

k-1 

implies 

2r 

^ < e for t^ <> t < » (2.4) 

k~l 

In the opposite case the state (2.1) is called unstable. 

The state is called asymptotically stable if it is stable and in addition 


2r 

lim Y 

t-*® L 

k=i 



(2.5) 
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This fundamental definition of stability by Liapunov has been refined and 
supplemented in various ways a n d reference should be made to the comprehensive 
texts by Minorsky [12], Krasovskii [13], LaSalle and Lefschets [14], and Hahn 
[15] 0 A host of "fine 11 definitions has been introduced, e 0 go, uniform stabi- 
lity, quasi-equiasymptotic stability, total stability, stability in the whole, 
etCo, [15] in order to classify possible effects of disturbances* 

For the purposes of the present monograph it appears to be sufficient to 
employ just three terms, namely, 

1) Asymptotic stability 

2) Marginal stability 

3) Instability 

Types 1) and 3) have been defined above<> Type 2) characterises a state which 
is stable, again as defined above, but not asymptotically stable* 

Expressed verbally, one can say that a state of equilibrium is asympto- 
tically stable if small disturbances, inflicted upon the system at a certain 
time, decrease with time* The state is marginally stable if the disturbances 
do neither decrease nor increase with time, and it is unstable if the dis- 
turbances increase with time* 

Side by side with Liapunov 1 s concept of stability, it is possible and 
meaningful to introduce alternate definitions * The two other most current 
ones are due to Poincar£ (orbital stability) and to Lagrange (boundedness of 
motions and orbits) , but the distinction vanishes for the special case of equi- 
librium states* Further, it would be of interest to examine the behavior of 
the system under continuous disturbances and under arbitrarily large disturb- 
ances,, For a discussion of the interrelation of various concepts of stability 
of dynamical systems, reference is made to Magiros [16]* Generalization of 
these concepts to continuous systems is not readily accomplished, because the 
notion of a metric has to be introduced, cf * Chapter V* 

Having accepted a definition of stability, the first step in the analysis 
of the state of equilibrium of a system involves the consideration of criteria 
which would permit to decide whether a given state is asymptotically stable, 
marginally stable or unstable * In dynamical, discrete systems two categories 
of criteria have been evolved, one being based on the construction of the so- 
called Liapunov's function (Liapunov's direct method), the other being based on 
the examination of solutions of equations of motion and, in continuous systems, 
modal expansions* In problems of stability of equilibrium the former is re- 
lated to the energy criterion which in turn, for certain systems, is equivalent 
to the static criterion (Euler method)* The latter is usually referred to as 
the kinetic criterion or the vibration criterion 0 For a detailed discussion of 
stability criteria reference is made to [5-7,14]* 

The applicability of stability criteria, as emphasized by Ziegler [1,5-7, 
17], strongly depends on the forces present in the mechanical system* If the 
forces depend explicitly on time, they are called instationary, if they do 
not, they are called stationary* The stationary forces generally depend on 
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both the generalized coordinates and generalized velocities * If velocity- 
dependent forces do no work in any elementary change of position, they are 
called gyroscopic forces (e.g* Coriolis forces); If they do negative work, 
they are referred to as dissipative (e.g* viscous damping, drag)* Among the 
velocity-independent forces, i*e*, forces which depend on generalized coordi- 
nates only, one encounters those which are derivable from a single-valued po- 
tential* These, such as for example, gravitational forces, are termed non- 
circulatory (or conservative)* All other velocity-independent forces are re- 
ferred to as circulatory, or nonconservative, or follower forces* Strictly 
speaking, dissipative, instationary and follower forces are all nonconservative 
forces, but the terms circulatory forces, follower forces and nonconservative 
forces are used in the literature with the same meaning and will be employed 
interchangeably In this work* 

The bulk of the present monograph is concerned with various classes of 
mechanical systems whose common feature it is that follower (or circulatory, 
or ,f nonconservative n ) forces are always present* To investigate the state of 
equilibrium of such systems, the analysis can be based on linearized equations 
of motion (or equilibrium, in certain cases) in the vicinity of the state to be 
characterized* Since follower forces are stationary, the system of equations 
obtained is autonomous (no explicit time dependence) and homogeneous (no forcing 
terms)* In discrete systems the circulatory nature of the follower forces mani- 
fests itself in that the force matrix is not symmetric, while in continuous 
systems the boundary value problem is nonself-adjoint* 

In this monograph both, extensions of Liapunov’s direct method and the 
"solution" method are employed, v?ith emphasis on the latter* In discrete sys- 
tems one is then concerned with a study of solutions of the type 



which leads to a study of the roots X^ (or oj^) of the associated characteristic 
equation [18]* If the real parts of all the Characteristic roots X^ are nega- 
tive (or the imaginary parts of all the characteristic roots are positive), 

the system Is asymptotically stable* By contrast, if at least one of the char- 
acteristic roots X^ is positive (or u)^ is negative), the system is unstable. 

If all the roots X^ are pure imaginary (or pure real), the system is margi- 
nally stable (critical case). Liapunov’s theorems assert that linearized ana- 
lysis is appropriate for asymptotically stable and unstable systems* In case 
of marginal stability of a linearized system, no statement can be made regarding 
the behavior of the actual nonlinear system* 

The nature of the roots X^ (or to^) can be determined without calculating 

the roots themselves* A variety of methods have been developed for this purpose 
[18], one of the most widely used being associated with the names of Routh and 
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Hurwitz.* Various definitions of follower forces as applied to continuous 
bodies have been discussed by Sewell [19], Nemat -Nasser [20 ], and Shieh and 
Masur [21]. 

For additional references, relating in particular the two areas of sta- 
bility and control, the reader is referred to the recent bibliography by Wang 
[108]. 


It has been called recently to the author's attention by P. C. Parks 
[106] that it was Hermite [107] who has established considerably earlier the 
conditions generally known as the "Routh-Hurwitz criterion. 
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CHAPTER III 


NONDISSIPATIVE SYSTEMS WITH TWO DEGREES OF FREEDOM 


3.1 An Illustrative Example 

To illustrate some characteristic types of dynamic behavior, in the 
vicinity of an equilibrium configuration, of a mechanical system subjected 
to a follower force, let us consider a double pendulum. Fig. 3.1, composed 
of two rigid bars of equal lengths l t which carry concentrated masses m^* 2m, 

m. The position of equilibrium cp^=* cp^ 1 * 0 Is to be investigated when the 

system is subjected to a force P acting along the bars (positive when com- 
pressive). For this purpose a perturbed configuration (cp^ 0, c p 2 ^ 0, but both 

small) is investigated. In this position elastic restoring moments ccp^ and 
c(cp 2 ~cp^) are Induced at the joints and the direction of the load P is specified 
to form an angle with respect to the direction of the bar in the equilib- 
rium position. 

This system has been investigated by Ziegler [5] for the special case 
a m 1, which may be termed the case of purely tangential loading. In Ref. 

[22] the full range - ® < <* < «> has been examined. As can be easily verified, 
the load P (and thus the system) is conservative only for a 0. The system 
may be considered a two-degree-of- freedom model of a continuous cantilever. 

The analysis, restricted to a linearized formulation, consists in the de- 
termination of the two natural frequencies of free vibration as a function of 
the loading. For sufficiently small loads both frequencies are real and the 
system is thus stable under an arbitrary small disturbance, exhibiting bounded 
harmonic oscillations. As the load is increased, instability may occur by 
either one frequency becoming zero (static buckling) at the critical loading 
and then in general purely imaginary, or the two frequencies becoming complex, 
having passed a common real value at the critical loading (marginal stability). 
The ensuing motion under a supercritical force in the first case is nonoscil- 
latory with the amplitude increasing exponentially (divergent motion), and the 
critical load can be determined statically by the Euler method. In the second 
case the ensuing motion is an oscillation with a definite period but with an 
exponentially increasing amplitude, and the critical load cannot be found by 
the Euler method because no associated adjacent equilibrium exists. The first 
case could be called n static instability" in view of the behavior at the criti- 
cal load, and the second "dynamic instability." In aeroelasticity , however, 
analogous phenomena have been termed "divergence" and "flutter," respectively, 
[23,24], and we propose to employ this terminology in the sequel. 
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Lagrange's equations in the form 


Tt (f^) - !£ ' Q k tt ‘ T ' v - k ' 1,2> O’ 1 ’ 

are used to establish the linear equations of motion, in which the kinetic 
energy T is 

T - \ mi 2 [3^ 2 + 29^2 + q.,, 2 ] (3.2) 

the potential energy V of the restoring moments is 

V - \ c [29^ - 2 9l 9 2 + 9 2 2 ] (3.3) 

and the generalized forces (due to applied loading) are 
Q x - Pi-foj - <*P 2 1 

(3.4) 

Q 2 * Pi[(l-or)9 2 ] 


These forms lead to the equations of motion 

Smi 2 ^ + m£ 2 $2 + (2c-P£)q>^ + (aPj£-c)cp 2 =» 0 

+ m£ 2 cp 2 “ cc Pi + [c - (l-Qf)P£] q> 2 = 0 
which, in turn, stipulating exponential solutions of the form 

9 k - \e imt (k - 1,2) 

yield the frequency equation 

4 2 n 

P o uj - p 2 <u + P 4 = 0 

where 

P Q * 2m 2 ^^ 

P 2 = me 2 [7c - 2(2-0-) Pi] 

P^ = c 2 - (l-or)[3cPX - (Pj() 2 ] 


(3.5) 


(3.6) 


(3.7) 


(3.8) 
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The four characteristic roots will occur in pairs, the positive and nega- 

2 

tive roots of the two values of <d obtainable directly from the frequency 

2 

equation. For a negative u) one root describes an exponentially divergent 
2 

motion; u) “0 corresponds to neutral equilibrium, the appearance of an ad- 
jacent equilibrium configuration (static buckling, divergence). A complex 
2 

value of to yields one root describing an oscillation with increasing ampli- 

2 

tude (flutter). The system is thus stable only as long as both values of u) 

2 

are real and positive. We are interested in the manner in which u) varies 
with P for different values of a. This is accomplished by investigating the 

2 

curves of P versus real values of <i> . 

Expanding the frequency equation we find that it is a general quadratic 
2 

equation in uu and P t 9 of the form 

A(iu 2 ) 2 + B(uo 2 P£) + C(P£) 2 + D(to 2 ) + E(Pje) + F - 0 (3.9) 

2 

where the discriminant, B - 4 AC, is 

4m 2 / 4 [2(l-a) + a 2 ] (3.10) 

2 

Since this expression is always positive, the frequency curves (P versus u) ; 

2 

P, <ju , real) are all of the hyperbolic type. 

Except for degenerate cases, which shall be noted, there are but two 
general types of hyperbolas, with regard to orientation in the real ur, P l 
plane, that may be encountered. These two types, qualitatively, are of 
"conjugate 1 ' orientations. 

In the first type, each of the two branches of the hyperbola yields a 
2 

single (real) value of u) for every load and the two values never coincide. 
Instability may occur only in the form of divergence or divergent motions. 

2 

In the second general type of hyperbola, the two values of uu , for any 

2 

load producing real values of w , lie on the same branch of the curve. For 
each branch there is one critical load at which the two values coincide. 

2 

Regardless of the behavior indicated by the real values of uo on the branches, 
these two critical loads always bracket a single, limited range of the load 

2 

"between" the two branches of the hyperbola, for Which the values of u) are 
complex and the free motions are of the flutter type. Since the system must 
be stable for sufficiently small loads, these critical loads must be of the 
same sign for any given value of <y. 
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The solution of the frequency equation is 


0 ). 


1.2 


2 7c - 2(2-or)Pjt T f4(Pi) 2 [2q-(y) + ty 2 ] - 4cPjt (8-g) + 41c 2 l 1/2 


4m I 


(3.11) 


from which P versus cjd^ ^ can plotted for any of. We determine, first, 

the critical loads corresponding to coincidence of frequencies (occurring 
in the second type of hyperbola) by setting the discriminant equal to zero 

2 

in the equation for ^ • This yields, for the critical loads, in nondi- 
mens ional form, the equation 


n m (8-tt) ± f (8-or) 2 - 4ir2q-g) + g 2 ]} 1/2 

c 2[2(l-tt) + a 2 ] 


( 3 . 12 ) 


Real values of these critical loads are associated with the second type 
of hyperbola, complex values with the first. We wish to determine the tran- 
sitional values of <*. Thus, setting this discriminant equal to zero. 


(8-a) 2 - 41 [2 (1-or) + a 2 ] - 0 


(3.13) 


yields the roots <* 


tr 


0.345, 1.305. 


Substituting this equation into that for ^ yields 

7c - 2 Pi (2 -qf) T [2Pj6(8-cv) 


K,2 


41c] 


(3.14) 


with o' = or 


tr 


tr 

0.345, 1.305. 


4m£ 


Thus two transitional values of a are obtained, at each of which the 
hyperbolas degenerate into two intersecting straight lines. Between these 
values of ot the second type of hyperbola is found to occur, and the phenome- 
non of flutter is thus limited to this range of <y. The corresponding critical 
loads are all positive (compressive). 

Consider next the constant term p, in the frequency equation. The Euler 

4 2 

method is equivalent to setting p^ =* 0(cc = 0), corresponding to intercepts of 
the hyperbolas on the P-axis. Setting p^ 3 0 we obtain for the Euler buckling 
loads, in nondimens ional form 




(3.15) 


12 



For real values of the load we must have or < | or tt > 1. Thus there are cri- 
tical values of <x(& cr ■ 1) , marking the limits of a range in which no ad- 

jacent equilibrium position occurs in the system for any value of the load. 


We note from the form of that the lower critical value of c* is a 

function of the elastic and geometric parameters of the system, and under a 
variation of these parameters might increase indefinitely, approaching one as 
a limit. Thus, for or < 1 there is a class of systems or loadings wherein the 
absence of an adjacent equilibrium configuration for any value of the load is 
a function of the elastic and geometric parameters. However, for a * 1, the 
terms in p^ involving P drop out entirely, leaving a positive definite ex- 
pression which contains the elastic parameter alone. Therefore, in this spe- 
cial case alone, we may say that it is the specification of the loading itself 
which results in the absence of any possible adjacent equilibrium configuration. 
We note here that in the case of a uniform continuous cantilever subjected to a 
load characterized by the same type of parameter, the Euler method reveals a 
similar critical value of the parameter, which in that case is one half. 


A third set of values of a of interest is denoted by &' > and is associated 
with a coincidence of an Euler load with a critical load for flutter. This 


occurs when a value of u) , at which id. 


is zero. Thus we set <d. 


0 , 


which is equivalent to setting p^ “ 0, p^ =* 0, simultaneously; i.e., 

[7c - 2(2-or )Pjfc] « 0 (3.16) 


and 


C 2 - (l-a)[3cPje - (P£) 2 ] * 0 (3.17) 

Solving the first equation for P£ and substituting into the second yields a 
quadratic equation in a, the roots of which are found to be a - 0.423, 1.182. 

In the sequel we restrict our detailed attention to 0 £ a < 1, as this 
range is somewhat more meaningful physically and is sufficient to demonstrate 
a connection between the various instability phenomena. Figure 3.2 shows the 
frequency curves for the various values of a of particular interest in this 
range. From these curves, in which both branches of the hyperbolas and their 
asymptotes are shown for completeness, we can determine by inspection the par- 
ticular character of the frequency curve for any a in the range 0 <> or £ 1, and 
we proceed now to a discussion of the behavior of the individual curves in this 
range and the various stability phenomena that they illustrate. 

For 0 £ of < c* tr the hyperbolas are all of the first type, and the behavior 

is as previously discussed. The frequency curves and the characteristic beha- 
vior of the nonconservative systems are qualitatively indistinguishable from 
the conservative case. Obviously, the Euler method would yield the lowest 
buckling load, which here marks the boundary between the single stable and 
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unstable ranges of the loading. A kinetic analysis would yield nothing addi- 
tional. With increasing values of a in this range, the hyperbolas draw closer 
to their asymptotes and finally degenerate into two straight lines at a » 

as previously noted. This case marks the first occurrence of a coincidence of 
the characteristic roots. 

For values of or greater than fl* tr in this range, hyperbolas of the second 

type* with the conjugate orientation, occur and pull away from their asymptotes 
with increasing c*. The upper branch lies entirely in the second quadrant, cor- 
responding to divergent motion, and such an instability follows flutter with 
increasing load. 

In this range, for or < or < ot* the coincidence of frequencies on the 

tr 2 

lower branch occurs at negative values of (jd , with divergent motion already 
characterizing both modes. Thus in this class of systems the boundary between 
the single stable and unstable ranges of the loading parameter is marked solely 
by the appearance of an adjacent equilibrium configuration in the first mode, 
and is obtainable by the equilibrium approach. The system is unstable for all 
higher loads. The critical loads corresponding to coincidence of frequencies 
do not mark any bound between stability and instability. 

Thus, for such systems, the Euler method would yield the critical load 
with regard to stability, even though the phenomenon of flutter is possible at 
some higher loadings. Conversely, the sole use of the kinetic method, if em- 
ployed so as to determine merely the critical loads corresponding to the coinci- 
dence of frequencies, would lead to erroneous conclusions. 

2 

For o' = ot' 9 to at the coincidence of the frequencies on the lower branch 
is zero. The sequence of instabilities with increasing load is the same as in 
the preceding range of a . 

For a > cf / the coincidence of frequencies occurs at positive values of 

and this critical point now marks the bound between a stable and unstable range 

of the load. However, for or ; < a < a the lower branch still intersects the 

cr 

load axis, and the two corresponding critical loads, both occurring in the first 
mode, now bracket a separate range of instability through divergent motion. 

Such a system is rather remarkable in that it displays, for different loads, 
losses of stability by both divergence and flutter. 

Thus for ex' < a < we have a rather interesting sequence of free motions 

with increasing load, resulting in multiple regions of stability and instability. 
This is illustrated in Fig. 3.3 by the frequency curve for the arbitrary value 
of a * 0.5. Such a system has characteristic free motions which include succes- 
sively stable oscillations, divergent motion, stable oscillations, flutter, and 
then divergent motion again for all higher loads. In such a situation the lowest 
critical load marking the appearance of an instability would still be a buckling 
load, obtainable by the Euler method. However, the existence of the second range 
of stability, above the second "buckling" load, as well as its upper limit, would 
be revealed only by a detailed kinetic analysis. 
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For a * or cr the two buckling loads, bracketing the lower region of insta- 
bility, coincide and the frequency curve is tangent to the load axis. Thus in 
this case there is a divergence instability at that isolated load, with no 
associated divergent motion for neighboring loads. The sequence of instabilities 

is otherwise the same as for a' < a < a . 

cr 

For ce ^ a? cr the lowest critical load was always a buckling load, obtainable 

by the equilibrium approach, and varied continuously with ot. Above a no ad- 

cr 

jacent equilibrium configurations occur and the lowest critical load is that 
corresponding to the coincidence of the frequencies. There is thus a disconti- 
nuity (jump) in the magnitude of the lowest critical load, at a ■ <* cr - Systems 

in the class a? cr < of £ 1 possess a single range of stability and of instability, 

with a sequence of characteristic free motions of stable oscillations , flutter, 
and finally divergence for all higher loads. 

The foregoing discussion could be extended to the values of a outside the 
range 0 < a £ 1 , but is omitted here for the sake of brevity. However, a plot 
showing the variation in all the critical loads for a wider range of a 9 including 
the entire region of flutter instability, and with the asymptotic behavior of 
the critical loads for divergence clearly indicated for extreme values of a, is 
given in Fig. 3.4. 

Considering systems corresponding to given values of a, this plot illus- 
trates that in systems displaying multiple regions (and types) of instability 
under compressive loading, the lowest critical load may correspond to either 
divergence or flutter. Also illustrated is the existence of systems displaying 
instability by divergence for both compressive and tensile loads. 

With the aid of the parameter <y in the simple model analyzed here, we have 
attempted to show a connection between instability phenomena of divergence and 
flutter by demonstrating a generic relationship between such disparate frequency 
curves as those characterizing a * 0 and <* * 1. Thus, such curves (and systems 
characterized by them) may be seen to be not of a singular or isolated nature, 
but part of a continuous "spectrum" of frequency curves. 

The justification for considering the entire range of - 00 < <* < + ® may be 
made clearer through the following observation. The type of loading specified 
may be considered as the result of a superposition of two component loads, cor- 
responding to constant-directional vertical loading (a“0) and tangential loading 
(a* 1), the two being kept in a constant ratio as the loading is varied. In 
such a perspective, 0 < ot < 1 corresponds to these component loads having the 
same sense. Then, ot < 0 and a > 1 corresponds to these component loads having 
opposite senses, with their relative magnitudes determined by the magnitude and 
sign of ot, and with positive load always corresponding to a resultant compressive 
loading. 

The effect of weights of the masses has not been included here, but our in- 
vestigations indicate that for small such constant loads the principal effect 
consists in shifting the frequency curves in the positive (negative) direction 
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of the abscissa for a suspended (inverted. Fig, 3.1) model. Referring to the 
frequency curve in Fig. 3.2 for a * 1, we can see that the shift caused by sta- 
bilizing constant forces would result in an intercept of the upper branch of 

2 

the hyperbola on the to * 0 coordinate axis. This is, in fact, the case ana- 
lyzed by Ziegler [5], in which the Euler method yielded a higher critical load 
than the kinetic method, and which has contributed to the general discrediting 
of the applicability of the static approach in nonconservative problems. This 
particular case is somewhat equivalent to the situation occurring herein for 
1 < a < c* tr » i n which, under compressive loading, the system becomes unstable 

through flutter, with the higher critical load, for divergence, of no conse- 
quence . 


3.2 General System with Two Degrees of Freedom 
3.2.1 Governing Equations 


Let us now generalize the specific results obtained in the previous section 
and consider a general system with two degrees of freedom. Let q^, q^ be the 

principal coordinates of the system and the equilibrium configuration q^= 0 

is to be investigated with regard to stability. The system is characterized by 
inertia (masses and m^) and by restoring spring constants and k^. Further, 

it is subjected to follower forces whose magnitude is dependent on a single pa- 
rameter. The linear equations of motion may be then written as 


+ a *2 P< l2 “ 0 
® 2 ‘4 2 + ^2^2 + a 21 Pq l + a 22 Pq 2 “ 0 


where < 2 *^ are assumed to be given. 

With the abbreviations 

*i 2 = k ±/ m ± 

“ij 

the equations of motion take on the form 


^1 + “i + a H Pq l + “l2 Pq 2 " 0 


q 2 + w 2 q 2 + a 21 Pq l + a 22 Pq 2 “ 0 


(3.18) 


(3.19) 


(3.20) 


We wish to characterize the position of equilibrium for various ranges of P and 
for various ranges of the system parameters. For this purpose we investigate 
solutions of the type 
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q k 


V 


i-u)t 


k - 1,2 


(3.21) 


which lead to the homogeneous set 

(-U) 2 +U) 1 2 4q' 11 P) A x + Qfi2 PA 2 “ 0 

“21 PA 1 + (" u,2+a) 2 2 ’ to 22 P ^ h * 0 

and finally to the characteristic equation 

* - S 2 - [(?i + ?2 ) + (o? ll' ta 22 )P ] 5 + a 12 a 21 p2 * ° 

where 

I * 0) 2 » $i - t»i 2 i - 1,2 


(3.22) 


(3.23) 

(3.24) 


This equation represent# a curve of second degree (a conic section) and may 
be written in the normal form 


F = a§ 2 + 2b|P + cP 2 + 2d? + 2eP + f = 0 


where the coefficients are given by 

a * 1; b - - I2-, c = a i\ a 22 ~ a 12 a 21 

d = ‘ (?i+ 5 2 >/ 2 ; • = («n5 2 - f «22 ? l )/2; f " ? 1?2 


(3.25) 


(3.26) 


The invariants of the characteristic curve are 


a b d 

bee 
d e f 


“l2 a 21 (? 2‘ 5 l ) /4 


(3.27) 


and 


a b 

b c 


- t(or 11 -a 22 ) + 4 o' 19 q' 91 ]/4 


12 21 

If A J* 0, equation (3.25) represents a regular second degree curve, namely, 


(3.28) 


an ellipse for 6 > 0 
a parabola for 6=0 
a hyperbola for 6 < 0 
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while, if A * 0, the curve degenerates into two real or imaginary straight 
lines. 


The system may lose stability, as we have seen, either by oscillating 
with increasing amplitudes (flutter) or by moving to another position of equi- 
librium (divergence). The critical values of P are associated with stationary 
points of the characteristic curve for £ > 0 (flutter) and with points of in- 
tersection with the P-axis (divergence). The ranges of system parameters 

and will be determined in which either flutter or divergence or both flutter 

and divergence may occur. 


3.2.2 Parameter Ranges 


A. 


a 12 a 21 


> 0 


It will be shown 
termine the points of 
hyperbola) with the P- 

that in this case no flutter can occur, 
intersection of the characteristic curve 
■axis which are 

First we de- 
(which is a 

p --('*V 

1 e Z -cf J /c 

for c f 0 

(3.29) 

P = - f/2e 


for c * 0 

(3.30) 

The discriminant 




D = e 2 - cf 

2 

3 ^ a ll ? 2‘ Qf 22 ? l ) 

+ ^12^21^ 1^2^^ 

(3.31) 


is positive and therefore there exists at least one real point of inter- 
section. 

To find stationary points of the characteristic curve F(§,P) =» 0, we 
have to calculate 

dF/dP = - (BF/d?)/(dF/dP) (3.32) 

and set this derivative equal to zero. If dF/dP ^ 0, it will be sufficient 
to examine 


dF/d? - 2§ + 2bP + 2d * 0 


(3.33) 


This equation is to be solved for § and a substitution made into the 
equation for the characteristic curve which in turn, solved for P, yields 

P =■ - [(e-bd) i^e-bd^-te-b^f-d 4 ") ] /(c-b 2 ) (3.34) 

In terms of system parameters the discriminant is 
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(3.35) 


(e-bd) 2 - (c-b 2 ) (f-d 2 ) - - « 12 « 21 (5 2~h )2,A 


Let us assume first that the two natural frequencies of the system are 

distinct. (The special case b* discussed separately.) This 

implies that the discriminant in Eq. (3.34) is negative and thus no real 
points with a horizontal tangent exist, indicating impossibility of occur- 
rence of flutter. 

For c i 0 and e < 0, the solutions P of the characteristic equation are 
positive. The system is marginally stable for all values of P in the range 
- <x> < P < P^, where is the smaller value of P from Eq. (3.29) and unstable 

for P £ P^ (see Fig. 3.5a). Similarly, for c > 0 and e > 0, both solutions 

of Eq. (3.29) are negative. Hence, the system is marginally stable for 
< P < ® (where |p^| < |P 2 |) and unstable for P £ p^ (see Fig. 3.5c). If 

c < 0, the two roots of Eq. (3.29) have different signs. In this case the 
region of marginal stability is given by P 2 < P < P^, while the system is un- 
stable for P £ P^ and for P £ P 2 (see Fig. 3.5b). 


If 0^2 and/or * 0, the characteristic equation takes on the form 


F s [or u P - (5-S 1 )3C* 22 P - (S-? 2 >] - 0 


(3.36) 


which represents two straight lines which may be considered as the limiting 
case of the hyperbolas of the previous section approaching their asymptotes. 
Again flutter cannot occur and stability can be lost by divergence only. The 
regions of marginal stability and instability are given in Fig. 3.6. 


In the special case * ***12 * 0 the eigenvalue curve degenerates Into 

two straight lines parallel to the P-axis. Thus no instability can occur for 
any value of P. 


C. 


“I2 a 21 


> 0 


a) Existence of Flutter 

In this subsection the ranges will be established in which flutter may 
occur or cannot occur. We solve Eq. (3.33) for P 

P - - (?+d)/b b + 0 (3.37) 

and substitute into the characteristic equation with the result 

A§ 2 - 2B§ + C - 0 (3.38) 
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which has the solutions 


5j TJ = (b ± J B 2 -AC ) /A for A * 0 

| = C/2B for A = 0 

where 

A =■ c - b 2 = - [< a 11 -° f 22* 2+ 4of 12 a 21^ 4 

B = Be - cd = - [(a' 11 -or 22 ) ( Q! 11 ? 2 “ Qr 22 ? 1 ) + 2a i2 a 21 (Sj+SplM 
C = cd 2 - 2bed + b 2 f 

and 

J B 2 -AC = (a 11 -haf 22 )(? 2 -l 1 ) J~ a i2 a 2l /4 

We first consider the case A ^ 0. With the notation 

2 

a ll " x ’ ff 22 * y ’ *12*21 * ~ p ’ v " § 2 /? l 

we obtain 

*l 9 H - B ± » - § 1 [ (x-y) (vx-y)-2p 2 (l*+-v)±(x+y) (l-v)p]/4 

which may conveniently be written in the form 
N I,II * ' ^1®1®2^ 4 

where 

gj^ 1,11 = y - x ± 2p 

g 2 1,11 = y - vx =F (v+l)p 

Similarly we write A in the form 
A = - hjh 2 /4 

where 

h x =* y - x - 2p 
h 2 = y - x + 2p 


(3.39) 

(3.40) 

(3.41) 

(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 
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Thus we find 


(3.49) 


- I, II I, II ,. 

5 I,II 5 1 8 1 80 /h ’ h 


12 


Since and are linear functions of x and y, it is easy to specify con- 
ditions for which 5^ or respectively, are positive or negative. The 

result is given in Fig. 3.7 in which the shaded area indicates that both 
solutions are negative and flutter cannot occur, while in the remaining 

area at least one 5^ is positive and flutter may occur. 

In the case A ® 0 we had § « c/2B which may be written as 


§ - g 2 r g 2 “ / [(x+y) (x-y) (v-1)] 


(3.50) 


The solid lines in Fig. 3.8 show the region where 5 > 0, while the broken 
lines indicate 5 < 0. 

It remains to investigate what happens if b * 0. The characteristic 
equations degenerate in this case to 

F e S 2 + cP 2 + 2d? + 2eP + f - 0 (3.51) 

and 

dF/d§ « 2§ + 2d =* 0 (3.52) 

It follows that 

5 =■ - d > 0 (3.53) 

and the equation for P is now 

cP 2 + 2eP + f - d 2 =■ 0 (3.54) 

whose solutions are 

P =■ - [e ±y e 2 -c( f-d 2 ) ] /c for c # 0 (3.55) 

P = - (f-d 2 )/2e for c ■ 0 (3.56) 

Since 

e 2 - c(f-d 2 ) = - a 12 « 21 (? 2 -? 1 ) 2 /4 > 0 (3.57) 

there is always a real solution for P which means that flutter may occur 
for b = 0. 
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b) Existence of Divergence 

The discriminant D of equation (3.29) may be written as 


D 



'lV 4 


(3.58) 


where 


* y - vx + 2p /v 
t 2 - y - vx - 2 P /v 


(3.59) 


It can be shown that t * 0 are equations of the tangents to the hyperbola 

2 1 
c 5 xy + p =» 0. 

Real points of intersection with the P-axis exist if D 2 0. The plain 
area in Fig. 3.9 indicates the region where D > 0, i.e., where divergence is 
possible, while the shaded area indicates that D < 0, where divergence is not 
possible. 

For the case c * 0 we have to use equation (3.30) and find that the 
solution P is finite except for the points T and T # at which the tangents 
t i * 0 touch the hyperbola c » 0 (see Fig. 3.9). Fig. 3.10 combines Figs. 

3.7 and 3.9 and shows a complete plot of the regions in which flutter or di- 
vergence, respectively, may occur or cannot occur. Fig. 3.10 can easily be 
constructed if p * “Q^o^l and v * ^2^1 are givetu 

3.2.3 Summary of Results 

If the loading parameter P is increased or decreased from its initial 
value (which need not be zero), the system may either remain stable, or it 
may lose stability by flutter or divergence. The results of the corresponding 
analysis are summarized in Fig. 3.11. Fig. 3.12 indicates qualitatively the 
ranges of stability for the loading parameter P for the various regions of the 
system parameters given in Fig. 3.11. Since substituting for or^ and 

-a^ for <*22 onl y reverses the orientation of the P-axis, the ranges of stabi- 
lity for the regions of system parameters indicated with a prime are obtained 
by substituting -P for P. Thus, Figs. 3.11 and 3.12 give a complete plot for 
the ranges of stability of the system for the case <*i2 a 21 < 

In many problems the initial value of the force P is zero, and one is 
only interested in how stability is lost first if a positive (or negative) 
force P is applied and increased. Fig. 3.13 shows whether stability is lost 
by flutter or divergence, respectively. 


22 



3.2.4 Special Case - u >2 

If the two natural frequencies of the system coincide, the characteristic 
curve degenerates into two real or imaginary straight lines. The characteris- 
tic equation (3.23) may then be written as 

F = [cP + b(b+ )(s-S 0 )][cP + (b- 7b 2 ^c)(5-§ o )] = 0 (3.60) 

for c # 0 

F s (2bP + § - 5 0 )(5 - 5 0 ) - 0 for c * 0 (3.61) 

where 

{„-«!- S 2 «•«> 

It is easily seen that for a i2 a 2l > ® and for <x l2 a 2l * ° t ^ ie resu * ts 8* ven 
In Figs. 3.5 and 3.6 hold if one sets =* § 2 * The results for a i2 a 2l < ® 

are illustrated in Fig. 3.14. If the system parameters fall into the regions 
II, III or IV, only divergence may occur. For system parameters in region I, 
flutter will occur for each nonvanishing value of P, while for system para- 
meters corresponding to point Q no instability will occur for any value of P. 
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CHAPTER IV 


DESTABILIZING EFFECTS 


4.1 Introduction 


It has been discovered by Ziegler [5] not quite two decades ago that in- 
ternal damping may have a destabilizing effect in a nonconservative elastic 
system. He considered a double pendulum with viscoelastic joints as a model 
of an elastic bar with internal damping and let a tangential force act at the 
free end. The critical loading obtained in complete absence of damping was 
found to be considerably higher than by Including damping at the outset of 
the analysis and then letting the damping coefficients approach zero (vanish- 
ing damping) in the expression for the critical force. 

This rather surprising and seemingly paradoxical finding was ascribed in 
later studies by Ziegler [6,7] to the possibility that internal damping is in- 
adequately represented by linear damping forces which are linear combinations 
of the generalized velocities and that the hysteresis effect should be taken 
into account. 

The destabilizing effect of damping was further elaborated upon by Bolotin 
[8] who considered a general two-degree-of -freedom system not related to any 
particular mechanical model and who found additionally that the destabilizing 
effect in the presence of slight and vanishing damping is highly dependent on 
the relative magnitude of damping coefficients in the two degrees freedom. 

Additional insight into the destabilizing effects of linear viscous 
damping in systems with follower forces may be gained by not merely applying 
stability criteria but by studying also the roots of the characteristic equa- 
tion (cf. Ref. [25]). Further, the results of the mathematical stability in- 
vestigations may be interpreted in physical terms by introducing the concept 
of degree of instability. It becomes then possible to carry out a gradual 
transition from the case of small damping to the case of vanishing damping 
and relate both of these cases to that of no damping. 


4.2 Illustrative Examples of Systems with Two Degrees of Freedom 
4.2.1 A Model 

For this purpose again a two-degree-of-freedom model is considered. Fig. 
4.1, composed of two rigid weightless bars of equal length £, which carry con- 
centrated masses * 2m, m^ * m. The generalized coordinates cp^, are again 

taken to be small. A load P applied at the free end is assumed to be acting 
at an angle (pure follower force). At the joints the restoring moments 

ecpi + b^ and c(q^-cp 1 ) + are induced. 
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The kinetic energy T, the dissipation function D, the potential energy V, 
and the generalized forces and are 

T “ \ ml 2 0* 2 + 24.^+ <(> 2 2 ) 

D - \ b^ 2 + \ b 2 (4» 1 2 - 2«p lt p 2 + cp 2 2 ) 

v - \ city 2 - 2q» lt p 2 + <p 2 2 ) (4.1) 

Q x ■ PjKcpj- cp 2 ) 

Q 2 - 0 


Lagrange's equations in the form 


d_ / 5X \ 3D_ _ dT_ 3V_ 

dt \a4» t / scp ± " ^cp ± acp t 


^i 


(i =■ 1,2) 


are employed to establish the linear equations of motion 

3mi 2 4> 1 + (b 1 +b 2 )tp 1 - (Pf-2c)cp 1 + wg 2 ^- b^ + (Pi-c)cp 2 = 0 

ml - b 2 <p^ - + mf $ 2 + b^ + ccp 2 = 0 

which, upon stipulating solutions of the form 

<Pi * A i e tt)t (i = 1,2) 

yield the characteristic equation 

P( n 4 + Pl n 3 + p 2 n 2 + P 3 n + p 4 = 0 

with the coefficients 

p o - 2 

Pj = Bj + 6B 2 

p 2 - 7 - 2F + B^ 

p 3 » B 1 + b 2 
P 4 = 1 


(4.2) 


(4.3) 


(4.4) 


(4.5) 


(4.6) 
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and the dimensionless quantities 


« - * <sr • 


*i- 


l (cm) 


1/2 


(i - 1,2) 


(4.7) 


F » — 
c 


In the absence of damping (®j“ B 2*0) , the characteristic equation Is a 
biquadratic 


2Q 4 + <7-2F)n 2 +1*0 


(4.8) 


4.2.2 Critical Loads 


From the assumed form of the time-dependence for the coordinates cp^ and 

on the basis of the kinetic stability criterion, it is evident that if all 
four roots of the characteristic equation are distinct, the necessary and 
sufficient conditions for stability are that the real roots and the real 
parts of the complex roots should be all negative or zero. In case of equal 
roots the general solution of cp^ will have terms which contain powers of t 

as a factor. If the real parts of equal roots are negative, the system will 
be stable (vibration with decreasing amplitude), but if these real parts are 
zero or positive, stability will not exist (vibration with increasing ampli- 
tude) . 

Let us turn our attention first to the case of an initially undamped 
system. The four roots of the biquadratic equation as a function of F are 
(a special case, at - 1, of the problem treated in Section 3.1) 

n i,2,3,4 - H* [' - (f ' /2 )] 1/2± [ F - (W 2 )] l/2 } <*■» 

which, depending on the values of F, may turn out to be pure imaginary roots, 
complex roots, or pure real roots. The nature of these four roots as F varies 
is graphically illustrated in Fig. 4.2 in which the values of the roots are 
given by the intersection points of the root curves and the horizontal plane 
which is perpendicular to the F-axis and passes through the given value of F. 
The illustrations in Fig. 4.2 include a perspective of the root curves, and 
also the orthographic projections on the real plane (Im O * 0), the imaginary 
plane (Re 0*0) and the complex plane (F * 0). 

It is found that there will always be two roots with positive real part 

if F > | - /2 3 2.086 =» F . For F = F there exist two pairs of equal roots 
^ e e 
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whose real parts are all zero. Thus the system is unstable for F > F fi . For 
F < F fi all roots are distinct and pure imaginary and thus the system is mar- 
ginally stable. 

We consider next a slightly damped system, assuming “ ^2 ** 0*01* No 

simple expressions for the four roots of the quartic equation exist; the 

numerical results obtained are illustrated in Fig. 4.3, where a perspective 

view is supplemented by three projections of the same three planes as in 

Fig. 4.2. Two roots will have a positive real part for F > 1.464 = F, . 

a 

Stability can be investigated directly without determining the roots of 
the characteristic equation by applying the Routh-Hurwitz criteria, which re- 
quire that all coefficients Pj(j = 0, ..., 4) of the characteristic equation 

and the quantity 

2 2 

x = PjP 2 P 3 • P 0 P 3 - Pi P 4 (4.10) 


be positive. For positive damping these stability conditions are satisfied, 
provided 


p 2 = 2[- F + | (7 +BjB 2 ) ] > 0 


X - 6B 2 2 ) {- F + M 


4B, 2 + 33BjB 2 + 4B 2 2 
2(B 1 2 +7B 1 B 2 +6B 2 2 ) 


(4.11) 


I B 1 B 2 ]} 


> 0 


For the system to be stable F must satisfy the following two inequalities, 
where g = 0 ^ g £ « : 

* < l+l B l»2 


(4.12) 


F < 48^33^4 + I 
2(p +7P+6) 


Since 


4B 2 + 33P + 4 < 7 
2 (8 2 +7p+6) 2 


(4.13) 


for whatever g in its range, it is evident that the critical load will be 
governed by the second inequality, i.e., 


F = 4B 2 + 33£ ± 4 1 

r J o T 9 D -1 D o 

2(p z +7p +6) * * 


(4. 14) 
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which depends on the ratio as well as the magnitudes of the damping coeffi- 
cients. 


For B^ « 1, as well as in the limit of vanishing damping, F^ becomes 


F 


d 


4B 2 + 336 + 4 
2(6 2 +76+6) 


(4.15) 


which is highly dependent on g and is in general smaller but never larger 

than F . The ratio of F, to F versus g is plotted in Fig. 4.4. It is noted 
e a e 

that when g * 4 + 5/2 =* 11.07, F,/F reaches its maximum value 1, The de- 

a e 

stabilizing effect is thus eliminated in this particular case, similar to that 
found by Bolotin [8]. For g * 0, F^/F e reaches its minimum value 0.16; i.e., 

the maximum destabilizing effect is about 84 percent in the present two- 
degree-of-freedom system. 

4.2.3 Case of Vanishing Damping 

The two disparate values of the critical load for no damping (B^ 33 0) and 
vanishing damping (B^-* 0) justify a more detailed investigation of the limiting 
process as the damping coefficients approach zero. 


Let us examine first the limiting process for the roots of the character- 
istic equation. It can be shown with the aid of the theory of equations [26] 
that if B^ « 1 and F < 4.914 this equation will have four complex roots. Let 

these roots be 


K ± iY 2 
l x i ± 


(4.16) 


Then one can write [26] 


2(Y l + 


P 

P 


l 

o 


4y i x i [(y i +x i )2+ (y 2 +\ 2 ) 2 ][(y 1 +x 1 ) 2 + (y 2 -x 2 ) 2 ] - J 3 

P o 

where p Q > P y , and X are as defined earlier. For vanishing damping 


(4.17) 


' Y 1 + X 1 = 0 

. Y i x i t<Y i +x i )2+ (y 2 +x 2 ) 2 ][(y 1 +x 1 ) 2 + (y 2 -x 2 ) 2 ] - 0 


(4. 18) 
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Hence 


or 


“V 


Y i 



- 0 


Thus 


Q 



± ix 2 


or 




f Y l ± Y 2 

l" Y l * iY 2 


(4. 19) 


(4.20) 


and a substitution of these four roots into the characteristic equation will 
show that they are the same as in the case of no damping. 

In the case of F ^ 4.914, the four roots will all be real for small B . . 
Let 1 


Q 




(4-21) 


In the limit of vanishing damping one can show similarly that either u^= v^= 0 
or u^= -v^ and u 2 = v^. For either alternative, substitution into the charac- 
teristic equation reveals that the roots are the same as in the case of no 
damping. 

Thus the conclusion is reached that whatever F the roots of the charac- 
teristic equation for no initial damping (B^= 0) are identical to those of 

vanishing damping (B^ 0). This implies that the motions of the system, for 

some given initial conditions, and whatever F, will be identical in the case 
of no damping (B^= 0) and vanishing damping (fi^ 0). 

We focus attention next on the loading F in the two cases and before 
passing to the limit consider small damping (B^ « 1). The positive real part 

of the roots of the characteristic equation in the range F^ < F < F g for seve- 
ral small values of B 2 and, as an example, B^ = 0 (i.e., 0=0) have been cal- 
culated and the results are displayed in Fig. 4.5, where F is plotted as a 
function of Re 0 for nine values of B^. This figure illustrates that for the 

larger values of B^, F^ represents the critical load because for F > F^ some 

roots will have a nonvanishing positive real part. A small increase of the 
load above F^ will result in a large increase of this real part. For small 

values of B 2 , however, even though F^ is still strictly speaking the critical 

load, its significance is lessened, because a small increase of the load above 
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F d will not result any longer in a large increase of Re Q. Large increase of 

Re Q will now be associated with small increase of a load which is slightly 

lower than F . For vanishing damping Re Cl * 0 for any F < F . We thus con- 
e e 

elude that during the limiting process the significance of F^ as a critical 

load is gradually transferred to F^, and at the limit of vanishing damping 

(B^ -* 0) F^ has to be considered as the critical load. It is apparent now 

that this conclusion could only be reached by considering the roots of the 
characteristic equation and not by merely applying the stability criteria of 
Routh-Hurwitz, Further, the reasons for the stability criteria yielding dif- 
ferent critical loads for no damping and for vanishing damping can be better 
understood by having considered small damping. 

4.2.4 Degree of Instability 

It was established in the preceding section that for vanishing damping 
(B^ -• 0) the four roots of the characteristic equation become identical to 

those of no damping (B^ =■ 0) while the stability criteria alone would in gen- 
eral yield disparate critical loads in these two cases. 

To establish a further connection between the mathematically derived 
critical loads for no damping (B^ =» 0) and vanishing damping (B^ -* 0) it ap- 
pears helpful to introduce into the discussion a concept which might be called 
"degree of instability" and which embodies a relaxation of the concept of in- 
stability as used when applying the kinetic stability criterion. According to 
this latter criterion a system is stable if a suitable disturbance results in 
a bounded motion in the vicinity of the equilibrium configuration; e.g. , the 
system is unstable if a disturbance leads to oscillations with increasing am- 
plitude (flutter instability). For this type of loss of stability one can state 
that from a practical point of view it will certainly matter how fast the am- 
plitudes increase. 

For example, should a suitable initial disturbance be merely doubled in a 
time interval which is large as compared to, say, some reference period, while 
the duration of the system being subjected to a nonconservative force is by 
comparison relatively short, the system may be considered practically stable, 
while, mathematically, of course, one would have to conclude that it is un- 
stable. 

In order to weaken the kinetic stability criterion, one could prescribe 
arbitrarily the allowable increase of the disturbance and would then obtain 
for a given value of the load a critical time, not unlike in the case of creep 
buckling. As an alternative, one could introduce another measure of the rate 
of amplitude increase. By analogy to decaying oscillations, where the loga- 
rithmic decrement serves the purpose of quantitatively assessing the rate of 
decay, we can use the same quantity also as a measure of the rate of amplitude 
increase. Thus 

A 

6 = log — S- (4.22) 

n+l 
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is 


where A is the amplitude of the oscillation at a certain time t and A , , 
n n+1 

the amplitude at t + T, where T is the period. In the present problem, neg- 
lecting the terms of decaying magnitude in the general solution of c p^, 6 gene- 
rally will be time- independent for flutter motions, except when the characte- 
ristic equation has equal pure imaginary roots. 

The kinetic stability criterion requires 6^0; i.e., A ^ A , . A neg- 

n n't" i 

ative 6 properly could be called the logarithmic increment and in a real system 
it is conceivable that 6 may attain a certain value & c in a certain interval 

of time without the system losing its stability in any practical sense. 

For 3 =* B^/B 2 “ 1 the critical load F is displayed as a function of 

B^ =* B 2 = B in Figs. 4.6 and 4.7. For however small but finite negative value 

of 6, the critical load for vanishing damping (B 0) will always be that for 
no damping (B = 0) , namely, F^ . However, the critical load for small damping 

(B < 1) may be smaller than F g but for finite 6, however small, is always 

larger than F d . For given 6 the value of (small) damping B which is associated 

with the minimum value of the critical load can be determined. 

For vanishing logarithmic increment (6 -+ 0) the function F(B) approaches 
a limiting curve which will contain the point F d on the ordinate. For 6 *» 0 

the stability region is closed; i.e., points on the curve 6 ** 0 in Fig. 4.7 
are stable, including the point F d on the ordinate. For B a 0 it is the point 

F fi which separates stability from instability, but belongs itself to the in- 
stability region. This limiting process provides thus additional insight into 
the generation of the critical load F d * 

4.2.5 A More General Model 


Further interesting types of behavior may be discussed if the follower 
force is generalized by means of the parameter oe as discussed in Section 3.1 
without damping. The system to be analyzed is that of Fig. 4.8 (cf. Ref. [27]). 
The kinetic energy T, the dissipation function D and the potential energy V 
are the same as in Section 4.2.1, while the generalized forces are the same 

as those in Section 3.1. The associated equations of motion are 


3m£, 2 $ 1 + (b 1 +b 2 )cp 1 - (PjJ-2c)9 l + - b 2 $ 2 + (aP2-c)9 2 =* 0 

(4.23) 

2 .. • 2 *« 

$i - b 2 4> 1 - ccp^ + mi cp 2 + b 2 9 2 - [(l-af)Pi - c] cp 2 =* 0 

which, upon stipulating solutions of the form (4.4) yield the characteristic 
equation 

P q C 2 4 + pjQ 3 + p 2 Cl 2 + p 3 0 + p^ - 0 (4.24) 
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with the coefficients 


Po- 2 


Pi - B 1 + 6B 2 


7+B. B. 


P 2 “ 2(2-0?) [-1+^] 

P 3 - d-«)(B 1+ 2B 2 ) [- P + (l-g( P+2 ) ] 


(4.25) 


1 / 2 . 


.1/2- 


P4 . (a . a){r .| [l+ (3g f )' ] } {F .| [l .(5Z !f) - ]} 

The dimensionless quantities Q, B., and F have been defined in Eqs. (4,7) and 
0 is again 0 = B ]/ B 2* X 

4.2.6 Root Domains of Characteristic Equation 

It was found that small damping rather than vanishing or large damping 
is the cause of the destabilizing effect, and thus only small damping (B « 1) 
will be considered in the sequel. 

Let us introduce first the following quantities: 

tj 1 12 

H " 6 Po p 2 " 16 p l 


~ 2 


3 (a-2) [f - 2 (2 -or) ] 


1 = P 0 P 4 ■ 4 P 1 P 3 + 12 p 2 

- 12 [4(a 2 -10a+10)F 2 + 4(25a-32)F + 73] 


J = 6 P o P 2 p 4 + 48 P 1 P 2 P 3 " 16 P o P 3 " 16 P 1 p 4 " 216 P 2 3 

= " 216 [( 8 “ 3+ 96a 2 -336cr+224)F 3 

- (348a 2 -1464a+1032)F 2 

- (1362a- 12 12)F - 161] 


(4.26) 

cont. 
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1 


K - p o 2 I - 12H 2 


s - 4[(a-l) 2 + 1] { F 


2[(of“l) Z +l] 


X {(8-of) + 6.325 [- (a-0.345)(a-1.305)] 1/2 }} 


{ 


X i F - 


2[(<*-ir+i] 


£(8-0f) - 6.325 [- (a-0.345) (a-1.305) ]) 1/2 |- 


A * I 3 - 27J 2 S | p 4 K 2 

X - Pl P 2 P 3 - P C P 3 2 - Pl 2 P 4 

- B 2 2 {(l-a)[8 2 + 128 + 4 - 8a(8+2) ] F 2 

- [2(8 2 +78+6) + 2(l-a)(3 2 +H8-10) 

+ (l-aXp^P+UJBjBj] F + [48 2 + 338 + 4 + (8 2 +78+6)B 1 B 2 ]} 
s b 2 2 £(l-a)[8 2 + 128 + 4 - 8a(8+2)] F 2 

- 2[8 2 + 78 + 6 + (l-a)(8 2 +118-10)] F + <4p 2 +338+4) } (4.26) 

where p o ,...,p^ and other symbols have been defined previously. 

It is known from the theory of equations [26] that: 

(a) When A < 0, the characteristic equation has two real and two complex 
roots. 

(b) When A > 0 and both H and K are negative, the four roots are all real. 

(c) When A > 0 and at least one of H and K is positive or zero, the four 
roots are all complex. 

These criteria lead to the different root domains shown in Fig. 4.9. 

The domain marked by crosses indicates the existence of four real roots; that 
marked by dots corresponds to two real and two complex roots; and that marked 
by horizontal dashes or by diagonal lines indicates the existence of four com- 
plex roots. The more detailed nature of the roots and the related stable and 
unstable behavior of the system may be deduced from the following. 

Domain A>0, H<0, K<0 

This domain is marked by crosses in Fig. 4.9. In it, p^, p^, and p^ are 
always positive; p^ is always negative. Applying the well-known Descartes' 
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rule of signs, regardless of the sign of p^, it is seen that in this domain 

the four real roots of the characteristic equation are always pairs of two 
positive and two negative ones. Consequently, this is throughout a region 
of instability by divergent motion. 

Domains A < 0 

These domains are marked by dots in Fig. 4.9. Let the two real and two 
complex roots in these domains be represented by 

f Pi * iP 2 

Q - ^ (4.2: 

l r l 1 r 2 

From the relations between roots and coefficients in the theory of equations 
[26] and the definition of the expression X in the Routh-Hurwitz criterion, 
the following relationships hold: 


2(p 1 +r 1 > - - \ (B^B^ < 0 

2[ Pl (ri 2 -r 2 2 ) + r l(Pl 2 + p 2 2 )] - - ^ - - \ P 3 

° (4.28) 

(p 1 2 + p 2 2 )(r 1 2 -r 2 2 ) - ^- |p 4 

4pi r i { [<.Pi*i> 2 + P 2 2 - r 2 2 ] 2 + 4 P 2 S 2 } " “*3 “ 8 X 

p o 

As p^ is always negative in these three domains, the third of the foregoing 
equations indicates that 

r/ > r x 2 (4.29) 


which, in turn, shows that the two real roots are of opposite sign. Hence 
these three domains are also regions of instability. Again, recalling that 
p^ < 0, it is seen from the foregoing four equations that the real part of 

the conjugate complex roots will be negative if X > 0 or if X < 0 and p^ < 0, 

but will be positive if X < 0 and p^ > 0. Whence it follows that divergent 

motion will prevail in this region, of the type as sketched in Fig. 4.10(a) 
if X > 0 or if X < 0 and p^ < 0, or as in Fig. 4.10(b) if X < 0 and p 3 > 0. 

It is noted that, if the system is undamped (B^ = 0), and r^ will vanish 

identically. The undamped system will therefore undergo divergent motion of 


35 



the type as sketched in Fig. 4.10(c). By definition, in all these cases, the 
system is unstable. 

Domain K > 0 


This domain is marked by horizontal dashes in Fig. 4.9. Let us denote 
the four complex roots in this domain by 


a 


Yi ± iv 2 
6 i 1 i6 2 


(4.30) 


Then, as before, the following relationships are obtained: 


2<Y!+V 


£i 

P o 


\ (B^WBp < 0 


(4.31) 


4y 1 6 1 [(y 1 +6 1 )2 + (v 2 +6 2 ) 2 ][(y 1 +6 1 ) 2 + (Y 2 -S 2 ) 2] = " 8 X 

p o 

which indicate that y ^ and 6^ (the real parts of the two pairs of conjugate 
complex roots) both will be negative if X > 0 but of opposite sign if X < 0. 


Now, within this domain, we have 
K S 8p 4 - p 2 2 > 0 


(4.32) 


or 

1 2 
p 4 >8 p 2 

which, in turn, leads to 


X < 


| (4p 3 - P lf2 ) 2 4 0 


(4.33) 


(4.34) 


Therefore, the real parts of the two pairs of conjugate complex roots are of 
opposite sign. The nature of these four roots indicates that in this domain 
the system will flutter. 

Domain A>0, H>0, K < 0 

This domain is marked by diagonal lines in Fig. 4.9. As the four roots 
are all complex, the signs of the real parts of the roots will also be governed 
by the signs of X as asserted in the foregoing. Thus the system will vibrate 
with decreasing amplitude (asymptotic stability) if the values of a and F are 
in those parts of this domain where X > 0. However, the system will flutter 
if the values of a- and F are in those parts where X < 0. 
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Further separation of stability from instability in the present domain 
is governed solely by the sign of X. This is illustrated for the four cases 
of 3 * 0, 1, 11.071, and », as shown in Figs. 4.11, 4.12, 4.13, and 4.14, where 
the regions shaded by diagonal lines are regions of stability; those shaded by 
horizontal dashes are regions of flutter; those shaded by small triangles are 
regions of divergent motion of the type shown in Fig. 4.10(a); those shaded by 
dots are regions of divergent motion of the type shown in Fig. 4.10(b); and 
those shaded by crosses are regions of divergent motion in which the time in- 
crease of the generalized coordinates is of the exponential type. 

It is to be noted that, in the present domain (A > 0, H > 0, and K < 0) , 
if the damping effects vanish, the four complex roots of the characteristic 
equation will all be pure imaginary and distinct. Thus the undamped system 
executes steady-state vibrations and is marginally stable throughout the domain, 
as found in [22]. 

4.2.7 Nature of Boundaries Separating Different Root Domains 

In this section, the boundaries given by X = 0, = 0, and K 3 0 will be 

examined. For the sake of convenience, the term "boundaries given by X * 0" 
will be restricted to mean only those parts of the curves given by X = 0 which 
lie in the domain A > 0, H > 0, and K < 0. 

Boundaries X(a,F,p) = 0 

On these boundaries, the characteristic equation has, by definition of X, 
two roots equal in magnitude but opposite in sign. These two roots are 

n i,2 < 4 - 35 > 


where p^ is positive for positive damping. It is found that the curves p^ = 0, 
p^ = 0, and X ® 0 have a common point of intersection which is given by 


20 2 + 50+2 

(4.36) 

F . F ' = + 1 
0 + 1 


Further, as p^ = 0 and X = 0 have only one point of intersection at (cx\ F ') 

on p^ =■ 0, it is evident that, along the boundaries given by X « 0, p^ is 

always positive. This can be seen from Figs. 4.11, 4.12, 4.13 and 4.14. 
Consequently, Q- 9 are two distinct pure imaginary roots. The sum of the 

other two conjugate complex roots is -p^/p^ e “ 2 P l* which is negative (for 

positive damping). Hence, along the boundaries given by X « 0, the charac- 
teristic equation has two pure imaginary roots equal in magnitude but opposite 
in sign and two conjugate complex roots with negative real part. Thus the 
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system will execute steady-state vibrations as a result of some initial dis- 
turbance. It is only in this case that the damped, nonconservative system can 
undergo such motions. 

Point of Intersection of X * 0, p^ * 0, p^ « 0 

At this common intersection point denoted by (#', F / ), the characteristic 
equation has two zero roots. The other two roots, being given by 

P o Q 2 + + P 2 - 0 (4.37) 

are two conjugate complex roots with negative real part. The two zero roots 
will induce two terms of the form c ^ + c^t in the general solution of cp i# 

Thus the system will execute divergent motion in which the increase of cp^ is 

linear with respect to time. This point (<*', F ') is the only one at which the 
stability region for the damped, nonconservative system is open. 

Points of Intersection of p^ ■ 0, X * 0, S «= 0 
Let us introduce the quantity 

S “ p l p 2 ‘ p o p 3 (4.38) 

such that 

X - p 3 S - Pl 2 p 4 (4.39) 

It can be shown that the curves p^ = 0, S * 0, and X * 0 have two points of 
common intersection, denoted by (cr*, F*) and (<**, F*"), where 



(5B 2 +228B+1440) ± £ (5B 2 +180B+800) 2 - 64006(8+6) 1 1/2 
16(158+112) 


F' 


5 (P+8 ) 

2(8+10-4a ) 


(4.40) 




_51§±8) 
2(8+10-4 a) 


These two points usually exist when £ is finite, but the point (a**, F //? ) ap- 
proaches infinity as £ -> ®. At the point (<**, F*), the characteristic equa- 
tion has one zero root, one positive real root equal to (“Po/p and two 

1/2 J 0 

negative real roots equal to "("^/Pq) and ~Pj/P 0 > therefore, the system 
will execute divergent motions. At the point (a", F'*), the four roots are 
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1/2 

one zero root, two pure imaginary roots equal to ± (-p^/pj) , and one nega- 

tive real root equal to -p^/p Q ; hence, after the Initial disturbance, the 

system will execute steady- state vibrations about a certain position which in 
general is not the position whose stability is being studied. 

Boundaries p^ * 0, Excluding Points (ot\ F 7 )* (cr # , F*), and (o'*, F*) 


Along these boundaries, the characteristic equation has one zero root and 
three other roots given by 

P 0 Q 3 + Pjfl 2 + P 2 ^ + P 3 ■ 0 (4.41) 

where, by the theory of equations and for small damping (B^ « 1), the three 

roots all will be real if p 2 < 0, but one real and two complex if p 2 2s 0. In 

the range of either F < F 7 or a > along p^ *= 0, the four roots are found 

to be one zero root, one negative real root, and two conjugate complex roots 
with negative real part. The nature of these four roots indicates that, after 
the initial disturbance, the system may execute transient vibrations and then 
come to rest at a position which in general is not the position whose stabi- 
lity is being studied. This phenomenon can be interpreted as a stabilizing 
effect of viscous damping because the same system with no damping would execute 
divergent motion. 

The curves p 2 “ 0 (i.e., H *= 0), p^ * 0, and K *= 0 have two common in- 
tersection points at (0.423, 2.219) and (1.182, 4.281). In the range of 
F 7 < F £ 2.2.9 along p^ = 0, the four roots are one zero root, one positive 

real root, and two conjugate complex roots with negative real part. In the 
range 2.219 < F < 3 along p^ = 0, the four roots are one zero root, one posi- 
tive real root, and two negative real roots. Thus, in the range of F 7 £ F < 3 
along p^ *= 0, the system will execute divergent motions. 

In the range F^<F £4.281 along p^ = 0, the four roots are one zero 

root, one negative real root, and two conjugate complex roots with positive 
real part; and thus flutter will occur. In the range F > 4.281 along p^ «= 0, 

the four roots are one zero root, one negative real root, and two positive 
real roots; hence the system will undergo divergent motions. 

~ 2 

Boundary K = 8p^ - p 2 “ 0 

The exact curve of K * 0 is 

K - 8P 4 - p/ - (p lPj - 1 PiV 64 "l 4 ) ■ 0 < 4 -‘ 2 > 
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As B i and, hence, and are assumed small, of the order of 10 , the last 

three terms in parentheses are higher-order terms and may be neglected * Thus 

K s Sp 4 - p 2 2 = 0 (4.43) 


is a boundary curve which is close to the exact curve K = 0* Substituting 
1 2 

g* p 2 for p^ in X, we have 

X “ " 8 ( p l p 2 _4p 3^ 2 * 0 (4.44) 


which indicates that the system will be unstable when a and F are on the 
boundary curve given by K = 8p^ - p^ = 0, except at the point where X vanishes 

and p- is positive (steady-state vibrations) * The instability mechanism, on 


the whole, will be of the flutter type, except at the points where the exact 
espressions of K and H are all negative (divergence) « 


4*2*8 Influence of Damning Ratio on Instability Mechanisms 

In the preceding sections, it was established that stability is possible 
only in the region (A > 0, H > 0, and K < 0), which is marked by diagonal 
lines in Fig* 4*9* In this region, the sign of X governs the type of motion, 
i*e., the system is stable if X £ 0 and unstable if X < 0. 


Critical loads for divergence, if any, are given by p, = 0; i*e. p they 

are 


div 


1 1 1 * e&o 


l/2j 


(4.45) 


On the other hand, critical loads for flutter, if any, are always given by 
X = 0, i.e., they are 


x 2(B 2 +9&-2)- tt(B 2 +llB-10)±(B+6)C(B 2 -22B+l)o 2 + 338a - 9S] 1/2 

flu 8(8+2) (a-1) (or-or ) 

o 

where 1 £ a r a , and 
o' 

„ 12 J3. + 4 

ff o 8(p+2) 


(4.46) 


(4.47) 


The two vertical lines a « 1 and <x c a Q (Figs* 4*11 to 4*14) are asymptotes 
of X « 0. For or *= 1, the critical load is given by 


F flu|o=l 


4B 2 + 33B + 4 
2(p 2 + 78 +6) 


(4.48) 
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which was studied in Sect, 4.2,2, For a - a 9 the critical load for flutter, 
if any, becomes 


4(8 + 2H4B 2 * 33B + 4) 

( f lu|o'=o' o p 4 + 7p 3 . 50p 2 _ 332p + 24 


(4.49) 


The curves of critical loads for p c 0 s 1, 11*071, and «> are illustrated in 
Figs* 4.11, 4*12, 4.13, and 4*14* 

For a « 0 (conservative case) in Fig* 4*11, the point (0, -1), which is 
an intersection point of two branches of the curves given by X c 0, is itself 
on the boundary given by X = 0; therefore, this point corresponds to steady- 
state vibrations of the system* The point (0, -1) is thus also a point re- 
presenting stability rather than a point which indicates an isolated critical 
load for the conservative system (or *= 0) with damping* However, depending on 
the ratio of damping coefficients, a nonconservative system (or £ 0) may have 
multiple critical loads for flutter, in addition to those for divergence, at 
the same value of a anywhere in the range or ^ 0, except for 5/9 < a ^ 1 where 
critical loads for flutter only will occur* Fig* 4*11 illustrates that, for 
P ® 0, flutter will occur for any or, except ct *= 0, while Fig* 4*12 shows that 
the smallest range of <x in which flutter is possible becomes minimum 
(5/9 < cl < 1*305) when the damping coefficients are identical (i*e*, p « 1)* 

It was found in Sect* 3*1 that the presence or absence of neighboring 
equilibrium positions was strongly influenced by the behavior of the noncon- 
servative loading and also by the constraints of the system* A further re- 
sult of this study is that the ratio of the damping coefficients may exert an 
analogous influence and may thus render the static criterion inapplicable for 
systems in which, without damping, the critical load could be determined sta- 
tically* For instance, it is seen that, in the range 1/2 < or < 5/9, the sta- 
tic stability criterion is applicable if p = » (see Fig* 4*14) but breaks 
down if p = 0 (see Fig* 4*11)* 

Similarly to applicability, the sufficiency of the static stability 
criterion (in the sense of supplying all critical loads) also depends on the 
ratio of damping coefficients* To exemplify this feature, let us examine 
again Figs* 4.11 and 4*14. In the range a < 1/2, we note that the static 
stability criterion is sufficient if p = « but proves to be insufficient if 
P « 0* The equation - 0 expresses, in fact, the static stability cri- 
terion, i*e., the condition of the static equilibrium of the system in the 
vicinity of its neutral configuration* Thus the static stability criterion 
is implied in the kinetic stability criterion, which is usually sufficient in 
determining all critical loads for the nonconservative system* 

It is possible to identify the range of a in which flutter cannot occur, 
and thus the application of the kinetic criterion is not required* However, 
this range will depend on the ratio of the damping coefficients. To determine 
this range, we consider the expression derived in this section* Flutter 

cannot occur if the quantity (p^- 22P + 1) c? + 33pa - 9p appearing under the 
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square root in that expression is negative* Thus flutter may occur in the 
following ranges: 


or 


where 


a ^ a 1 and o' £ if P > or p < 


<x l > a > ar 2 if a 1 > P > 


a l,2 " 11 ±/12 ° = 


121.954 

0.046 


and 


„ = -33 ft. ± IP. ( 3 ^ + 2?_7j +-36)Jl! - , a A* A a) 

*1,2 2(p-a 1 )(p-a 2 ) * U l ? P ? V 


(4.50) 

(4.51) 

(4.52) 

(4.53) 


If P - a^ or P = the range in which the kinetic stability criterion must 

be considered will be only or ^ 3/11 0 Consequently, if there exists any range 
of or which is outside the foregoing specified ranges, the static stability 
criterion alone will be sufficient to determine all the critical loads, de- 
spite the nonconservativeness of the loading* However, according to the pre- 
ceding section, if or < or 7 or a > o the static stability criterion defini- 
tely will be applicable but not necessarily sufficient in determining all 
critical loads* 


4.2*9 Possibility of Elimination of Destabilizing Effects 

Critical loads for flutter in the undamped system analyzed in Sect. 3.1 
are given by the equation K(or, F, B^) = 0 with the terms due to small damping 

neglected; i.e., by the equation 

K(or,F) = - [4 (a 2 -2ert-2)F 2 + 4(a-8)F + 41] = 0 (4.54) 

Critical loads for flutter in the damped system analyzed here are given by 

X(a,F,p) = 0 (4.55) 

whose loci constitute, in fact, a family of curves in the or - F plane with p 
as the parametric constant. Different curves of the critical load for flutter 
will be obtained if different values are assigned to P in X(<*,F,P) » 0. 

To study the interrelation between the curves of critical loads given by 
K(a,F) = 0 and X(cr,F,p) = 0, let us examine the envelope of the family of 
curves defined by X(or,F,p) = 0. It is known that, if an envelope exists, it 
must satisfy 

X(or,F,p) = 0 (4.56) 
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and 


|p X(or,F,p) - 0 (4.57) 

Elimination of P in these two equations yields 

(F-2)[(l-a)F-2][4(l-a)F-5] 2 *K(a,F) - 0 (4.58) 

where K(a,F) is as defined before* However, this equation may contain some 
curves which are other than the envelope* Deleting these, the true envelope 
is found as given by 

[(l-a)F - 2] • K(<*,F) - 0 (4.59) 

Thus the curve for critical flutter loads of the system with no damping is a 
branch of the envelope of the family of curves of the critical flutter loads 
of the same system with damping* This remarkable relation shows a significant 
connection between the two governing equations of the critical loads for 
flutter of the undamped and the damped systems. 


In consequence of the foregoing relation, it appears possible to elimi- 
nate the destabilizing effect of damping on the critical loads for flutter in 
the damped system if we choose the value of P which defines a curve of the 
family X(<*,F,P) * 0 tangent to K(a,F) ■ 0 (the envelope) at the given value 
of or* Eliminating F in X(or,F,p) ■ 0 and (B/dp)X(ar,F,p) * 0, we find that 
this value of p is given by the positive, real root of the quintic 


(8a - 3) (7a - 3) (4a - 3)p 5 - (896a 4 - 5,936a 3 + 8,196a 2 

- 3,870a + 594) p 4 - (12,800a 4 - 60,928a 3 + 82,680a 2 

- 38,664a + 5832)p 3 - (80,128a 4 - 365,280a 3 + 502,416a 2 

- 234,576a + 34,992)p 2 - (353,280a 4 - 1,480,320a 3 
+ 1,925,856a 2 - 874,800a + 128,304)P - (838,656a 4 

- 2,941,056a 3 + 3,411,072a 2 - 1,469,664a + 209,952) - 0 (4.60) 

and the critical load for flutter in this case is given by 

F . (15-32q)p 2 + (24- 128a) p + (84-496a) (4>61) 

2[(6-17a+8a 2 ) p 2 + (24-92af32a 2 )P + (120-484cH-256a 2 )] 


which will be identical to the critical loads for flutter of the same system 
with no damping. 

For example, if the elimination of the destabilizing effect of damping 
for the case or - 1 is desired, p must be equal to the positive, real root of 
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the quintic 


p 5 + 6p 4 - 86p 3 - 884p 2 - 2612p - 2448 - 0 

(4.62) 

i.e.. 


P = 4 + 5 /2 = 11.071 

(4.63) 

which, together with a ■ 1, yields 


F = “ - /2 ■ 2.086 

(4.64) 


The critical load for or =* 1 in the undamped system determined in [5,22,25] 
is identical to the value we obtained in the foregoing. The complete elimi- 
nation of the destabilizing effect for this case is thus attained, as is il- 
lustrated in Fig, 4.13, For a * 3/4, a similar procedure will show that the 
destabilizing effect is completely removed when p * ». This is illustrated in 
Fig. 4.14. 

The possibility of a complete elimination of the destabilizing effect de- 
pends on the existence of a positive, real root in the foregoing quintic. The 
range of a where the elimination of the destabilizing effect is of interest to 
us is, of course, 0,423 £ ot £ 1.305. However, it is found that in the range 

3/7 < a < 3/4 (4.65) 


the quintic has no positive, real root. Thus, in this range, the system will 
always experience some destabilization for whatever value of p in its range 
0 £ p £ 

For instance, let us consider the case a = 0.6, where the critical load 
for the system with no damping is 


F. =|g O 7 - 6 /5) = 2.033 


(4.66) 


while the critical load for the system with damping is given by 

_ . 1.4B 2 + 11.46 + 2 - (B+6)ro.36B 2 + 2.88B + 0.36l 1/2 

d (3.2|3 + 6.4) (a Q - 0.6) 


where 

= B 2 + 12 B + 4 
“o 8(p + 2) 


(4.67) 


(4.68) 


The ratio of F^ to F^ versus p is plotted in Fig. 4.15. It is noted that the 

value of F^/F^ increases as p increases and approaches 29/5(37-6 /5) ■ 0.984, 

instead of 1, as the upper limit when p approaches infinity; i.e., the desta- 
bilizing effect of damping is at least 1.6 percent if the value of a is kept 
at 0.6. 
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In the range 1.182 < a < 1.305, the undamped system has multiple critical 
loads for flutter given by K(cr,F) * 0. However, an investigation of the roots 
of the quintic shows that, for any ct in the range 1.182 £ a £ 1.285, there is 
only one positive, real root which defines a curve of the family X(a,F,p) * 0 
tangent to the lower part of K(a,F) * 0. Thus, in the range 1.182 £ ol £ 1.285, 
the damped system has no critical load which is given by the upper part of 
K(a,F) - 0. 

As an alternative, the possibility of eliminating the effects of damping 
could also be studied by equating the frequencies first and then the critical 
forces, obtained with and without damping. The frequency of the undamped 
system is given by 

Im n - | [7 - 2(2-a)F] 1/2 (4.69) 


while the frequency of the system with damping is given by 
Xm n - (i2) 1/2 r (W - (l-tt) (B 1 +2B 2 )F ,1/2 

' p i' L B i + 6b 2 J 


(4.70) 


Equating the two expressions and eliminating F in K(or,F) = 0 leads to 

28 (a - p 2 + 4(16a 2 - 33a + 9)p + 4(182a 2 - 297a +81) =0 (4.71) 

which, in turn, gives the range of Eq. (4.65) in which elimination of the 
damping effect is not possible for positive damping. 

Fig. 4.16 illustrates the function (3(c0 which insures elimination of 
damping effects. For completeness, the required values of negative (3 in 
the range 3/7 < ex < 3/4 have also been indicated. 


4 . 3 D amping and Gyroscopic Forces in Systems with Two Degrees of Freedom 


The joint effects of follower forces, linear viscous damping, and gyro- 
scopic forces (i.e., velocity-dependent forces which do no work) have been 
studied in Ref. [28]. Considered was the system with two degrees of freedom 


$1 + a nqi + a 12 q 2 + - 0 

«2 + a 21 q l + a 22*2 + b 21*l + b 22*2 " 0 


(4.72) 


The matrices a_ and b^ 
symmetrical part: 


can be resolved uniquely into a symmetrical and anti- 



(4.73) 
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where 


where 


*12 " *21 “ ( * 12 + a 21 )/2 > 


P ” ^*12'*21 )/2 


(4.74) 


and 


11 “l2l I “11 °12 


21 22 


n h , 


22 


r 0 o)"' 


+ < 


- 0 ) 0 


(4.75) 


b U = b 21 " < b 12 +b 21> /2 ‘ 


u> - ( b 12 - b 2l )/2 


(4.76) 


By a suitable transformation of the form 
**1 ” q l cos ^ ' ^2 8in 

(4.77) 

q 2 * sin <p - q 2 cos 

it is possible to make either a* 2 or b* 2 to vanish. Choosing the first pos- 
sibility and writing again q^ q 2> b 12 ... for q^, q 2 , b* 2> the following 
system of equations is obtained: 

^1 + *ll q l + pq 2 + b ll q l + (b 12 + “ 0 


q 2 - pq x + a 22 q 2 + (b 12‘ ® )q l + b 22 q 2 = 0 


(4.78) 


The system has a potential energy function (is noncirculatory) if p - 0, it 
is purely circulatory for a^= a 22 = it * s nongyroscopic for u) - 0, and is 

undamped if °* 


Solutions are sought in the form 

A Xt 
q i " A i e 

which lead to the characteristic equation 


c X 4 + c_X^ + c_X^ + c~X + c. =* 0 
ol234 


where 


c =1 
o 

C 1 “ b ll + b 22 


(4.79) 


(4.80) 


(4.81) 

cont. 
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(4.81) 


C 2 ’ a ll + a 22 + (b ll b 12- b 12 2) + " 2 


C 3 " a ll b 22 + a 22 b ll + 2ptu 
C 4 “ a ll a 22 + p2 


For stability it is required that £ 0 (i * 1,2, 3, 4) and that in addition 


X - CjC 2 c 3 - c o c 3 2 - C]L 2 c 4 > 0 


(4.82) 


If cy* 0 this additional condition takes on the form 


c 9 “ *c c y > 0 
2 o 4 


(4.83) 


It is to be noted that if 0 but cjt 0 (or 0 but 0) the first in- 
equality cannot be satisfied and thus the equilibrium is always unstable re- 
gardless of the actual values of the nonvanishing c^. 

Let us now examine the special case of an undamped system, b^** The 
coefficients are then 


c l « 0 , 


C 2 " a ll + a 22 + “ * 


c 3 = 2puj, c 4 


a ll a 22 + P 


(4.84) 


Since c^ a 0 and 0 the system is unstable regardless how small the 
follower (circulatory) forces and the gyroscopic forces are. 

Another special case of interest in which explicit results can be ob- 
tained is 


a ll " a 22 = a > °' 


b U " b 22 " b > 0, 


b !2 "° 


(4.85) 


Then we have 


= ]* 2b, c 2 = 2a + b Z + ou , c 3 = 2 (ab-fpco) , c^= a Z + p 2 


(4.86) 


For stability we must require 


ab + pu) > 0 

2 2 
ab 4- pojb - p >0 


(c 3 > 0) 

(X > 0) 


(4.87) 


It is again seen from the second inequality that no stability is possible for 
b = 0 or for small b. The damping coefficient b has to be sufficiently large, 
namely 
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b > pw + Jp 2 (t» 2 +4a) ^ / (2a) 


(4.88) 


to insure stability. In the absence of purely gyroscopic forces, a) = 0, the 
stability condition is 

b > p//a (4.89) 


4.4 Discrete Systems with Many Degrees of Freedom 

Generalizing the findings concerning destabilizing effects found with 
specific examples of systems with two degrees of freedom, it is possible to 
state a number of theorems which are applicable to a rather broad class of 
systems with N degrees of freedom, (Ref. [29]). In particular, it can be 
shown that not only slight viscous damping, but all sufficiently small 
velocity-dependent forces may induce a destabilizing effect. 

The system considered is assumed to be holonomic and autonomous, and 
is subjected to a set of generalized forces, = Q^(F); j =» 1,2,...,N, which 

are defined as linear functions of a real, finite parameter F. This para- 
meter, (0 ^ F < «►) , is associated with the magnitude of the externally applied 
forces, Qj = 0 for F = 0. 

Let 

q j * Qj ” 0, (j = 1,2,. ...N ; 

be the equilibrium state of the system. With M 
matrix, and 

N 

v ■ 1 1 Va <4 - ,l) 

J,k-1 


ij-dt) (4 - 90) 

- the generalized mass 


the strain energy function, assumed to be positive definite, the equations 
of motion of the undamped system may be written as 

M jk 9 k + K jk q k « Qj j,k = 1,2,...,N (4.92) 

where the summation convention on all repeated indices is implied and will 
be employed in the sequel. 


Let us assume that the generalized forces, Q ^ , are given as linear 
functions of the generalized coordinates 


^ * F V“k 


J.k - 1,2,. ...N 


(4.93) 


where K = [K . ] is a nonsyrame trie matrix, and F a real, finite parameter. 
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For F = 0, (4.92) represent the equations of free oscillation of the undamped 
system which we assume to possess N distinct, non-zero frequencies. 

In conjunction with (4.92) we shall consider the following linear system 
M jk q k + eG jk q k + *jk q k = Q j j (4.94) 

where e is an infinitesimal quantity, G =* [G^l a generally non-symmetric 

matrix with prescribed constant elements. For e * 0, Eqs. (4.94) reduce to 
Eqs . (4.92). 

In the following sections we shall prove that the critical load of system 
(4.92) is an upper bound for the critical load of system (4.94) when 0(e 2 ) can 
be neglected in comparison with 0(e). Only the effect of velocity-dependent 
forces on the critical load of the system for flutter will be considered. The 
effect of these forces on the critical load for divergence is discussed in 
Refs. [6,7]. 

In the present context, therefore, the theorems proved in the sequel are 
applicable only when a linear system loses stability by flutter. 

It is also of importance to note that an autonomous, linear, dynamic 
system can lose stability by flutter if and only if a solution of the form 

q^ = k = 1,2,...N, admits, at least, one od with negative imaginary 

part. Further, we will employ the well-known property of linear autonomous 
dynamic systems of the type (4.92) that the roots of the characteristic equa- 
tion are either real or pairs of complex conjugate numbers. 

Let us first consider the effect of slight viscous damping. Thus we as- 
sume that G * [Gj^] is a symmetric, non-negative matrix. 

We take solutions of (4.92) and (4.94) in the form q^ = i = /-l, 

and obtain 


- • 2 " j k*k + <V'V*k ■ 0 <4 - 55) 

- + eicuGj^ + (Kjk _ FKj k ) = 0 j,k = 1,2 N (4.96) 


Systems (4.95) and (4.96) are each a set of linear, homogeneous equations 
in A^. They have, therefore, nontrivial solutions if and only if the deter- 
minant of the coefficients of A^, in each set, is equal to zero. These con- 
ditions yield 

det | a l = 0 (4.97) 


det l a jk + eiuC jk! = ° 


(4. 98) 
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2 i 

vhere - - to + (Kj^-FK^), and det | a j k 1 denotes the determinant of 
the matrix [a^]. 

For F ■» 0, Eq. (4.97) yields the natural frequencies of the free vib- 
ration of the undamped system. We assume that these frequencies 


(u>. < t» 0 


2 . 
a) N ) 


are distinct and non-zero. We now increase F and assume that for a certain 
value of F, say F , Eq. (4.97) yields, at least, a double non-zero frequency. 
e 2 2 

Let us suppose that, for F =* F^, is equal to cd^ (see Fig. 4.17(a)), while 

all other (N-2) frequencies of the system are distinct and non-zero. If F is 
now increased beyond this critical value F^ , Eq. (4.97) will yield a pair of 

complex conjugate roots and, consequently, the system will oscillate with an 
exponentially increasing amplitude (flutter). We shall refer to F^ again as 

the critical load for the system without damping. 

Let us now consider Eq. (4.98). For F = 0, the roots of this equation 
are all located on the left-hand side of the imaginary axis in the complex 
ioo plane. As we increase F, at least, one of these roots approaches the im- 
aginary axis, and for a certain value of F, say F^ , Eq. (4.98) yields, at 

least, a real value for uo (see Fig. 4.17(b)). If F is now increased beyond 

this critical value F, , at least, one of the roots of (4.98) becomes complex 

a 

with negative imaginary part. The system, therefore, loses stability by 
flutter. We shall refer to F^ as the critical load for the system with 
damping. 


In the sequel we will first study a system with two degrees of freedom 
and then extend our results to more general systems. 

We expand the frequency equation of the damped system as follows 

det |a 1k + eiuC jk | = det |a 1k | + ciuX3 1k a kj - e 2 u> 2 det |G^ k | ; j,k - 1,2, (4.99) 


Jk 

kj 


jk 


Jk 


jk' 


where a 

assume that det |G. 


Is the cofactor of the element a^ in the det iSjjJ* Moreover, we 
^ 0 (the case of det 1 G 1 =0 will be discussed later). 


Then, for det 


k ! ' v v '“- jk * 

finite and e of infinitesimal order, we may neglect the 


last term on the right-hand side of Eq. (4.99) and obtain 

0 j,k = 1,2 


det I a jk1 + * iu} G jk akj 


(4.100) 


Theorem 1 . The critical load, F , is an upper bound for the critical 
2 6 

load, F^ , when 0(c) can be neglected in comparison with 0(c). 
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Proof . For F * F^, Eq. (4*100) has, at least, one real root, u) - to, and 

the other roots are either real or complex with positive imaginary parts (the 
possibility of complex root with negative imaginary part is excluded, as it 
contradicts the assumption that F^ is the critical load). Therefore, for 

F * F^ and au real, det | a jjJ and are bot ^ and we must have 

det l a jkl “ 0 (4.101) 


G Jk a 


kj 


j»k - 1,2 


(4. 102) 


However, det a 


F d S F e* 


jk' 


0 cannot admit real roots if F > F 


Therefore 


Let us note that F^ can equal F^ if and only if the real root of (4.102) 

can be made equal to the double root of (4.101) for F 3 F g . This, of course, 

depends on the other parameters of the system and may not always be achieved, 
as exemplified in Sect. 4.2. 

We now retract to Eq. (4.99) and consider the case when det |g,,| “0. 

JK 

The frequency equation of the damped system with two degrees of freedom is 
now given by Eq. (4.100), independently of the order of magnitude of c. 
Following the line of reasoning similar to that used in the proof of Theo- 
rem 1, we conclude that the critical load of the system without damping is 
an upper bound for that of the system with damping, no matter what the order 
of magnitude of c may be. Therefore, we state the following theorem. 


Theorem 2 . The critical load, F^, of the system without damping with 
two degrees of freedom is an upper bound for the critical load, F^, of the 
system with damping for all finite values of c when det |g , | 3 0. 

The proof of Theorem 1 was an immediate consequence of a property of 
the frequency equation of the system with damping and with two degrees of 
freedom. The problem becomes more complicated if the system has more than 
two degrees of freedom. However, one may still use a similar line of rea- 
soning. 

We expand Eq. (4.98), collect the terms of like power in e, and obtain 

det i a j k + G j k l 3 |a 1 + eiu) Gj fc a k ^+ 0(c 2 ) + ..., j,k-l,2,..N (4.103) 

The first term on the right-hand side of this equation is a polynomial of 
2 

degree N in u) and may be written as 

det | a jk | s p(u) 2 ) s P n o) 2N + P^u) 2 (N-1) + ... + P q (4.104) 
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k 1 2 

Similarly, the term G^a » which is a polynomial of degree (N-l) in 00 , can 

be written as 


V“ 3 * “to 2 ' ■ «H-i 


tu 


2 (N-l) 


V: 


,00 


2 (N-2) 


+ ... + R 


(4.105) 


Therefore, Eq. (4.103) becomes 

det | a j k + eioo G^ k | = P(cd 2 ) + iew R(uj 2 ) + 0(c 2 ) + . 
2 

We neglect 0(e ) and higher in Eq. (4.106) and obtain 


(4. 106) 


P(u> 2 ) + lew R(tu 2 ) = 0 


(4.107) 


for the frequency equation of system (4.94). We now set 00 = X + iey and 
2 2 

substitute into P(u> ) and R(cu ) to obtain 


P (X+ieY ) = P(X 2 ) + lev f 2X d -%- *1 + 0(e 2 ) + ... 

L d(\ ) J 


(4. 108) 


R (X+ieY ) - R(X ) + 0(c) + ... 
Therefore, Eq. (4.107) becomes 

. 2 , 


P(X 2 ) + icX 2y dP ^ - ) + R(X 2 )1 + 0(c 2 ) + 
L d(X^) J 


= 0 


(4.109) 


Neglecting terms of order higher than e, we must have 


P(X 2 ) =0, Y = - , P'(X 2 ) = dP ^ ' * 0 

2P'(X‘ £ ) d(X ) 


(4.110) 


The first equation in (4.110) is the frequency equation of the system 
without damping and the second equation defines, to the first order of appro- 
ximation in e, the effect of slight damping on the frequencies of the system. 

2 

The constraint given by P ; (X ) ^ 0 indicates that the perturbation method 

2 

breaks down when P(X ) 53 0 admits double roots. For F = 0, the roots of equa- 
2 

tion P(X ) = 0 are all real and distinct. Thus in this case, to the first 
order of approximation in e, the roots of Eq. (4.107) are 


\ = ± \ + ieY 


k* 


R < x k z > 

2P'(X k 2 ) 


1,2,. ..,N 


(4.111) 


2 2 

where X^ are real, positive roots of P(X ) * 0. The system can perform 
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oscillations only with an exponentially decaying amplitude and, therefore, 
all k ° 1,2,...,N are positive, real numbers. 


We shall now assume that the system with damping is stable for all 
F < and consider the following cases: 


(a) F < F, < F 

N ' d e 

(b) F e < F < F d ; F e < F d 


( 4 . 112 ) 


2 

For case (a), P(X ) * 0 yields N distinct roots. From Eq. (4.111) we 
then obtain k ** 1,2,...,N, which are, by our assumption, all positive, 

real numbers. 


2 

For case (b) , P(X) * 0 has, at least, one pair of complex conjugate 
roots. We denote these roots by -** *S) anci from (4.111) obtain 


+ 


0) 


1,2 


± (a ± ip) - ic 


*r(« ± ip) 2 i 

2P'[(a ± ip) 2 ] 


(4.113) 


which indicates that, for F > F^, the system with damping admits, at least, 

one complex frequency with negative imaginary part. This, therefore, con- 
tradicts the assumption that the system is stable for F^ > F fi . We are thus 

forced to take F, £ F in order to remove the contradiction, 
a e 

Let us note that, for F « F d * F^, Eqs. (4.111) can be used only for the 

distinct roots of P(X^) = 0. The perturbation method, which was introduced 

2 2 

here breaks down if P'(X ) “ 0(c) while R(X ) is non-zero. We shall not, 
however, concern ourselves with a detailed study of this case here and simply 
admit the possibility of F^ “ F e * l n fact, as > F e renders the system 

unstable, we can only conclude that F^ ^ F e * Therefore, we may state the 
following theorem. 


Theorem 3 . The critical load, F^, of system (4.92) is an upper bound 
for the critical load, F^, of the system with slight damping when c is suffi- 
ciently small. 

For an arbitrary specified matrix G “ [G^] due to any type of velocity- 

dependent forces (including gyroscopic forces), system (4.94) may become self- 
exciting. That is, for an infinitely small value of F, the frequency equation 
of this system may possess complex roots with negative imaginary parts. In 

these cases we shall agree to define F, “ 0 as the critical load of this system. 

d 

On the other hand, the frequency equation of system (4.94) may yield roots 
with only positive imaginary parts for F » 0(c). This indicates that this 
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system is stable for small values of the load parameter F. However, as we 
increase F, one of these roots may move toward the imaginary axis in the iu) 
plane. Therefore, for a certain value of F, say F d , the frequency equation 

of system (4.94) may yield a non-zero, real root. In this case, if we then 
increase F beyond this critical value F^, the frequency equation will have a 

root with negative imaginary part and the system will flutter. We shall re- 
fer to as the critical load of system (4.94). On the basis of the above 

preliminaries it is now possible to follow the same chain of arguments out- 
lined previously and establish the following more general theorems. 

Theorem 4 . The critical load, F^, of system (4.92) is an upper bound 
for the critical load, F^, of system (4.94) when c is sufficiently small. 

G * need not be a symmetric, positive definite matrix. 

Theorem 5 . The critical load, F fi , of system (4.92) with N ■ 2 is an 
upper bound for the critical load, F^, of system (4.94) for all finite values 
of e when det | G | * 0- G * [G^] need not be a symmetric, positive definite 
matrix. 

From the above results we immediately conclude that, in a linear system 
with N degrees of freedom, subjected to nonconservative (i.e. circulatory) 
forces, not only slight viscous damping but all sufficiently small velocity- 
dependent forces have, in general, a destabilizing effect. Moreover, the 
critical load, F^, is highly dependent upon the structure of the matrix 

G ° [G^] but is always bounded from above by the critical load F^. This in- 
dicates that, even at the limit as c -» 0, F^ is in general less than F fi . Let 
us explore this point in more detail for a system with two degrees of freedom. 

For c finite, the steady state motion of the system is possible if the 
frequency of the oscillation satisfies the following equations (see Eq. (4.99)) 

det |a jk | - eV det |6 jfc | - 0, G jk a kj - 0; j,k - 1,2 (4.114) 

In this case, one may solve the second equation in (4.114) for to as a function 

of F and then substitute the result into the first equation to obtain a re- 

lationship between F and c. In this manner a stability curve, in the F-e 
plane, may be constructed (see Fig. 4.18). However, from Theorem 1 we im- 
mediately conclude that, in general, thJ } curve suffers a finite discontinuity 
at c ■ 0, This means that, although for e *= 0 the critical load is F , for 

4 - ^ 

c “ 0 the critical load is given by which is, in general, less than F^. 

Therefore, the point F^ is, in general, an isolated point in the F-* plane 

(Fig. 4.18). This phenomenon was interpreted physically in Sect. 4.2. 
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4.5 Dest abilizing Ef fects in Continuous Systems 


4.5.1 Introduction 

It was shown in Sect. 4.4 that In a general circulatory system with N 
degrees of freedom not only slight viscous damping, but all sufficiently small 
velocity-dependent forces, such as Coriolis forces in vibrating pipes conveying 
fluid, or other gyroscopic forces, may have a destabilizing effect. 

For a continuous system, however, which possesses an infinite number of 
degrees of freedom, no such theorems are as yet established. To study the ef- 
fect of viscous damping forces in such systems most investigators, in general, 
reduce first the continuous system to a discrete one by means of, for example, 
the Galerkin method, and then study the reduced, discrete system [8,30,31]. 

But, as was shown in Sect. 4.4, a discrete system does, in fact, always have 
this property, except in very particular cases. Therefore, by this approach 
one does not know whether the original continuous system also exhibits the 
same behavior or whether it is produced only through the reduction procedure. 

Let us show that the presence of sufficiently small velocity-dependent 
forces in a continuous elastic system subjected to follower forces does, in- 
deed, have a destabilizing effect (cf. Ref. [32]). To this end, a cantile- 
vered, continuous pipe conveying fluid at a constant velocity is considered. 

The internal and external viscous damping forces are also included, and then 
it is proved that the critical flutter load of the system may be reduced by 
almost 507. for some combinations of these velocity-dependent forces. The 
method of analysis effectively reduces a complicated nonself-ad joint boundary 
value problem (without discretization) to a simple frequency analysis by uti- 
lizing fully the fact that the velocity-dependent forces are sufficiently 
small. 

It is of obvious interest to test the accuracy of the widely used Galerkin 
method with a two- term approximation. It is to be noted that such an analysis 
of this approximate method, for the case when the equations of motion of the 
system also contain mixed time and space derivatives, has been carried out in 
Ref. [32] for the first time. 

Critical flutter loads of the system, for small velocity-dependent forces, 
and also for large values of Coriolis forces, were obtained by using the Galer- 
kin method with a two- terra approximation. The results are then compared with 
the exact solution. It is then shown that the two- term approximation yields 
sufficiently accurate values for the critical flutter load only if the velocity- 
dependent forces are small. Thus, for large values of Coriolis forces the cri- 
tical load obtained by the Galerkin method with a two- term approximation may be 
greatly in error. 

4.5.2 Cantilevered Pipe Conveying Fluid 

We consider a cantilevered, uniform pipe of length L and internal cross- 
sectional area A, conveying fluid at a constant velocity U. A nozzle whose 
opening is n times smaller than A is placed at the free end of the system, as 
is shown in Fig. 4.19. 
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We shall assume that the material of the pipe obeys a stress-strain re- 
lationship of the Kelvin type, i.e., 

a = Ee + Tie (4. 115) 


where E is the modulus of elasticity and T| is the coefficient of viscosity. 
Under the assumption of plane sections remaining plane, the moment -curvature 
relationship, for small deformations, is 


m_ = . ra!z + D a 3 y 1 
EI Ux 2 E a x 2 at J 


(4.116) 


where M is the resultant moment at section x and at time t, I the moment of 
inertia, and y the transverse deflection of the pipe. With u denoting the 
displacement in the x direction, and z the distance of each fiber from the 
neutral axis, we also have 


Hz 
I 9 


du 

dx 



The equation of motion may now be stated as 


(4.117) 



(4.118) 


where p is the resultant lateral force exerted on the pipe. This lateral 
force may be decomposed into three parts. The first part is due to the 

S 2 y 

inertia forces and is given by + (m + m- ) — £ , where m is the mass of the 

1 dt 

pipe per unit of length, and m^ the mass of the fluid contained within the 
pipe. The second part is due to Coriolis acceleration and is given by 

a 2 v 

+ 2ra^U d x d t » and finally, the third part, which is due to equivalent com- 
pressive force induced by the flux of momentum out of the pipe, and is given 
2 d 2 y 

by + m n U n — £ . Therefore, the equation of motion becomes 
1 dx' 2 


*2„ *2 .2 
^ = (m4m ) ^ + 2m, U ^ 
dx z 1 dt 


2 

r dxdt + m i u n ri 

ox 


(4.119) 


and substitution from (4.115), (4.116), and (4.117) into (4.119) finally 
yields 

ei zi + m rfe - + m i u2n rl + 2 m i u & + rf - 0 < 4 - 120 > 


dx 


dx dt 


dx 


dt 
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If we include also the effect of external damping in the form K , 
is a constant, and introduce the following dimensionless quantities: 


where K 


5-J ■ 


2 2 
m.U nlT 

- F Z 


El 


V El ’ 


m + m. 


P' . 


E(m+m^)L 


k 2 l 4 


El (m+ra^ 


C4.121) 


then we obtain 

A. A 


^ + 6' + F 2 2_| + 2 / F |=|- + Y' P- + H » 0 

2>S 4 a ? 4 a T a ? 2 v n d|dT 3t a T 2 


(4. 122) 


To study the effect of small viscous damping forces and Coriolis forces, 
we now let 


6' - v6. 


2vy , and 




(4.123) 


where v is a small parameter. The equation of motion, (4.122), and the boun- 
dary conditions at § ■ 0, 1, may then be written as 






ar a T 


a§ aT 




at § » 0 


(4.124) 


0; .t { - 1 

H a§ 

We wish to study the stability of system (4.124) when \> is sufficiently small. 

We let y * ^r(5)e^ uyr , and reduce (4.124) to the following boundary value 
problem 

+ F 2 — - u) 2 + + iuw j^6 + 2pF ^| + 2vt^ = 0 


* = * ' = 0; at 5 = 0 


(4.125) 


0; at ? 
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where prime denotes differentiation with respect to 5. 

We then set f » e^; X * X + iva, and obtain 

(X+iva)^ + F^(X+iva)^ - 0)^ + ivti) [6(X+iva)^ + 2£F(X+iva) + 2y] “ 0 


(4.126) 


which is the characteristic equation of system (4.125). Expanding (4.126) in 
a series of powers of v, we are led to 

{x 4 + F 2 x 2 - (I) 2 } + (iv) {4aX 3 + 2F 2 Xa + id(6 X 4 + 2pFX + 2y)j- + 

+ (iv) 2 |6a 2 X 2 + F 2 a 2 + u)(46aX 3 + 2pFa)j- + (iv) 3 -[4Xa 3 + 66<ua 2 X 2 ]- + 

+ (iv) 4 ja 4 + 4u» 3 xJ- + (iv) 5 (u» 4 ) « 0 (4.127) 

2 

Next, we equate terms of like powers in v, neglecting 0 (v ) and higher, 
and finally arrive at 


2,2 




cu 


(4.128) 


U) 


6X + 28FX + 2y 


X - X + iva 


The 


2\(2X + F) 

4 

V X *5 

solution to system (4.125) may now be written as *(£)=■ ) A^e J , where 


J-l 


; j * 1,2, 3, 4, are constants which can be obtained from the boundary con- 
ditions at 5 =* 0,1. That is, they must satisfy the following four linear, 
homogeneous equations: 

iv 


ih 

j=i 


Aj = 0 


(4.129) 


J-l 

y x. 3 a.. Xj -o 

L j j 

J-l 
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System (4*129) has non-trivial solutions if and only if the determinant 
of the coefficients is identically zero, i.e., the frequency equation is of 
the form (written out explicitly in [32]) 

A « A (X.) - 0 (4.130) 

J 

This relation may be rewritten with the aid of (4.128) as follows, after 

2 

expanding it in terms of powers of v, and neglecting 0 (v ), 

A * { F^ + 2u > 2 + 2 w 2 ch cos + F 2 <u sh ^ sin X^ } - 

" lv { ( 2 ^2^ 2 )[ (X l 4 ‘ X 3 4) + (X 3 3 " 3X 1 V ch X 1 sin X 3 + 

+ (3\ 1 X 3 2 - X x 3 ) sh cos X 3 ] + 

+ 6 X^X 3 ^ ch Xj^ cos Xg - ^X^Xj 2 sh Xj^ sin X 3 + 2\^\^ sh Xj^ sin X 3 + 

+ 2 X 1 3 X 3 3 sh X x cos + Xj 2 ^ 4 ch X x sin X 3 - Xj ^ 2 ch X x sin X 3 ] + 

4 

+ ( SX 3 + 2y SX^/) + 6 X x 3 X 3 2 ch X x cos X 3 + 

2 X 3 ( 2 X 1 2 -tf 2 ) 

+ 4Xj 2 X 3 3 sh Xj^ sin Xg - 2 x^X 3 sh X x sin X 3 - 2X^7^ ch Xj_ sin ^3 + 


/ ^X^ + ^Y \r 5 4 

+ 511 ^ + 


+ X 1 2 X 3 4 sh Xj cos X3 - X^X^j 2 sh Xj cos X^ j- - 0 


2 F 


♦vg-w. 


(4.131) 

(4.132) 


The first term in braces, in Eq. (4.131) is the frequency equation when 
v * 0 , and the second term, to the first order of approximation in v, indi- 
cates the effect of small viscous damping forces and Coriolis forces. For 
v ** 0 , we obtain the frequency equation of a purely elastic cantilevered beam 
subjected to a compressive force which stays tangent to the axis at the free 
end. The critical value of the load, in this case, is F 2 * 20.05, which was 
first computed by Beck [33]. 

For non-zero but sufficiently small values of v and for small F, all the 
roots of equation (4.131) are located to the left of the imaginary axis in 
the complex ico plane. As we increase F, at least one of these roots approaches 
the imaginary axis, and for a certain value of F, say F^, Eq. (4.131) yields 

one purely imaginary root itu ■ iu> c . If we now increase F beyond this critical 
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value F^, one of the roots of (4.131) becomes complex with negative imaginary 

part, and the system oscillates with an exponentially increasing amplitude. 
Therefore, for given values of 6, 0 and y, we shall seek critical values of 
cu » (« c (real), and F * F d which identically satisfy (4.131). This is illus- 
trated in Fig. 4.20 where, for 6*1, 0*1, and y * 0, real (A,) tnd imagi- 

2 1 

nary (-Z^) P arts of A are plotted against the values of a) . Similar results 
may be obtained for other values of 6, 0, and y. 

It may also be of interest to establish the destabilizing effect of 
Coriolis forces, internal viscous damping forces, and external viscous damping 
forces independently. 2 

F d Z 

To this end, we let 6 * y * 0, 0*1, and with y d * — obtain, from Eq. 

TT 

(4.131), = 1.78. Similarly, for 0 * y = 0 and 6*1, the critical load is 

obtained to be y^ * 1.107. However, for 0=6=0 and y = 1 we get y d = 2.035, 

which is equal to the critical load of the system when no velocity -dependent 
forces are present. That is, although sufficiently small Coriolis forces and 
internal viscous damping forces have a destabilizing effect in this continuous 
system, external viscous damping forces do not have the same effect. 


The combined effect of velocity-dependent forces on the value of the cri- 

p 2 

F d 

tical parameter y d = — is shown in Figs. 4.21 and 4.22. In these figures 


TT 

the parameter y^ is plotted against 0/6 for various values of y. The horizon- 
tal dashed line in these figures represents the critical value of y^ when no 

velocity-dependent forces exist and the cantilevered column is subjected to 
a compressive follower force at the free end (Beck's problem [33]). 


It is important to note that the stability curves shown in Figs. 4.21 

and 4.22 have a finite discontinuity at v = 0. That is, although for v = 0 

2 2 + 2 
we have F = F^ = 20.05, for v = 0 , the critical value of F is, in general, 

less than 20.05. 


It may also be of interest to explore the order of magnitude of v for 
which the destabilizing effect of velocity-dependent forces still exists. 

This may be accomplished by considering v large and seeking values of cu and 
F for which Eq. (4.130) is identically satisfied. We note that, in Eq. (4.130), 

Vjj j * 1,2, 3,4, are defined as functions of w and the other parameters of the 

system through Eq. (4*126). In order to circumvent the difficulty of solving 

polynomials with complex coefficients, we let 6 * y * 0 and put X = iT| in Eqs. 
(4.126) and (4.130). 
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I 


The critical values of ou and F may now be evaluated by a computer. The 
computer may be instructed to obtain the roots of Eq. (4.126) for given para- 
meters, and then calculate A, (Eq. (4.130)). These results are shown in Fig. 

« 2 

F d / a ' 

4.23, where y d * — is plotted against values of J — 9 by a solid line. 

TT 

The dashed line in this figure corresponds to the critical when the Galerkin 
method with a two-term approximation is employed for the analysis as follows. 


We consider a set of orthonormal eigenfunctions, {cp n (§)}> obtained by sol- 
ving the following eigenvalue problem n 






A" " * n <Pn 


(4.133) 




dc P n 

dg 


0 ; 


at ? 


(4.134) 


at § 


(4.135) 


We 


.2 A 3 

d 'Pn d % 

di 2 V 

oo 

then let y * ^ ^n^ T ^n^* substitute it into the first equation in 
n*l 


“ * r-i 

(4.124), multiply both sides of this equation by 6y ■ ) cp m ( D ^^( t) » and 

m*l 

integrate the result from zero to 1 with respect to 5 to obtain 


— r 5 + Y (to 2 6 + F 2 b ) a + 

,2 L m mn mn ^m 

dT , 

m=l 

00 

+ v I < 6 % 26 mn + 2f5a mn + 2y W * °* n = 
m 

where 

tp m - cosh x a | - cos X^g - a m (sinh X^g - sin X^g) 

sinh X^ - sin X Jn 
°m cosh Xjg + cos X^ 


(4.136) 


(4.137) 
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a 

an 


(4.137) 


r 1 dtp -» « 

■ J d r *n d * 

o 

N . 2 “ »« 


b * 
Bin J 

,1 d 2 cp 

l J* 


I 

4 

1 ; for 

6 - 
mn J 

1 W 15 * 

o 

0 ; for 


System (4.136) is a set of nonself-adjoint, linear, second order, homo- 
geneous, ordinary differential equations which admit solutions of the form 

« A m e ituT . To obtain the critical values of F^, we seek conditions under 

which ou becomes complex with negative imaginary part. System (4.136), how- 
ever, consists of infinite number of equations each with infinite number of 
terms. This, therefore, leads to a determinant which possesses an infinite 
number of rows and columns. 

It is quite common to let m,n * 1,2 in Eqs. (4.136) and reduce this sys- 
tem to only two linear, homogeneous differential equations [8]. Hence, the 
characteristic equation becomes a polynomial of degree four, which can easily 

2 

be solved. The values of F , which render at least one real root and all the 
other roots complex with positive imaginary parts, are then taken to be ap- 
proximations to the critical flutter loads. 

For sufficiently small values of v, we may neglect terms associated with 

2 

v in the characteristic equation, and using Routh-Hurwitz criteria, calculate 

2 2 

approximate values of the critical load F * y^rr . In Table I these approxi- 
mate flutter loads are compared with the exact values obtained in the previous 
section. From this table we observe that, for sufficiently small v, the 
Galerkin method with a two- term approximation yields very accurate results. 

2 

We note also that, for v * 0, this approximate method gives F « 20.15 as com- 

2 

pared with the exact critical load, F * 20.05. 

The above conclusion, however, does not imply that, for v finite, the ap- 
proximate method should necessarily give sufficiently accurate results. In 
fact, as is shown in Fig. 4.23 for 6=7=0, the critical flutter load ob- 
tained by the approximate method (dashed line in Fig. 4.23) can be greatly in 
error for relatively large values of the Coriolis forces. We note that, for 
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^ smaller than 0.25, the resulting error, when the Galerkin method with 

a two-term approximation is used, is less than 5 percent and decreases as the 
value of v decreases. 

Among the other studies concerned with the destabilizing effects of 
velocity-dependent forces (and in particular linear viscous damping), mention 
should be made here of the papers by Leipholz [34] and Leonov and Zorii [35]. 
In Ref. [36] Bolotin and Zhinzher have used an expansion in fractional powers 
of the damping parameters and have established the conditions under which li- 
near viscous damping has no effect on the critical load for flutter. By con- 
trast, Zorii [37] was interested in determining the maximum effect which 
(small) linear viscous damping may have on the critical load. 


4.6 De stabilizing Effects Due to Phenomena Other than Linear Viscosity 

4.6.1 Thermoelastic and Hvsteretic D amping; 

Not only linear viscous damping, but other types of dissipation mecha- 
nisms are associated with destabilizing effects. In Ref. [38] a general for- 
mulation of the stability analysis of elastic continuous systems subjected to 
follower forces in the presence of thermomechanical coupling was presented and 
applied to the problem of a cantilever under a tangential follower force at 
the free end. A pronounced destabilizing effect of thermoelastic dissipation 
was found to exist. Bilinear hysteretic damping was studied in Ref. [39] where 
it was shown that it may have a destabilizing effect similar to linear viscous 
damping, but that this effect disappears for a large class of hysteretic sys- 
tems. 


4.6.2 Magnetic Damping in a Discrete System 

Damping in a system can be realized also through the interaction of a 
current carrying conductor with a magnetic field. Leibowitz and Ackerberg 
[40] have found that the motion of an electrically conducting, perfectly 
flexible wire placed in a transverse magnetic field will also be damped, but 
in a manner somewhat weaker than the familiar viscous damping. 

It is of interest to examine the effect of such magnetic damping on the 
stability of equilibrium of some circulatory elastic systems, cf. Ref. [41], 
where additional details are given. A simple system with two degrees of free- 
dom is considered first, and a destabilization is found to be caused by the 
magnetic field. 

The system consists of two rigid weightless rods, each of length A, car- 
rying concentrated masses m and 2m and acted upon by a follower force P (Fig. 
4.24). The rods OA and AB constitute portions of electrical circuits having 
resistances R^ and R^, respectively, and are constrained to undergo at most 

plane motion. A uniform magnetic field (jj^ acts in a direction perpendicular 
to the plane of possible motion. 


63 



A displacement from the equilibrium configuration (cp^ ° ^ will re- 

sult in elastic restoring moments ccp^ and c (q> 2 " <Pj) at t ^ ie hinges, and motion 

of the system in the magnetic field will induce a potential difference between 
any two points of the rods given by 

JjC x * d t (4.138) 

where the integration is taken over the conducting path joining the points and 
v c v(s)is the velocity of the conductor. The potential difference will re- 
sult^in the generation of a current, according to 




(4.139) 


and therefore a force per unit length of conductor given by 


f 

rj 


e x 


© 


o 


(4.140) 


where e is a unit vector in the direction of the current# The force distri- 
bution"^. 140) will of course be normal to the conductor and in a direction 
which opposes the motion. 


For the system being considered the distributions f^ and f^ are 


f i “ r i?i 

f 2 - r 2 ( 29 x + q> 2 ) 


(4.141) 


where 

r. -® o 2 i 2 /2R. (j * 1,2) (4.142) 

and the dots indicate differentiation with respect to time t. Taking as ge- 
neralized coordinates the (small) angles cp^ and cp^, the kinetic energy T and 

the generalized forces Q^, are found to be 
T =* mjfc 2 (3q> 1 2 + 2 ^ 92 + /2 


Q x - (PX-2c) 9l - (Pi-c)cp 2 - (r L + 4r 2 )X 2 ^/2 (4.143) 

Q 2 * c ( 9 1 "9 2 ) “ r 2 ^ 2 ^ 2 ^l + ^2^ 2 

These quantities are substituted into Lagrange's equations to obtain the 
linear equations of motion: 
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3m£ 2 cp x + ( r i + ^ r 2 ^ 2 ^l/ 2 “ ( p 4-2c)cp^ 

+ mJi 2, cp 2 + + (Pi-c)cp 2 *= 0 

2 m 2 • 2 

mi $]_ + r 2^ “ cc Pi + $2 

+ q>2^ 2 + Cc ?2 " 0 

The general solution of the system (4.144) is taken in the form 

<P j - X A j k ^ (j " 1 » 2 > 

k=l 

and leads to the characteristic equation 
4 3 2 

q o <« + q^o) + q 2 ci) + q3U) + q^ ° 0 

with the coefficients being given by 

o 2 fl 4 7 2 
a = 2m £ /c 
o 

q x “ (r^+ 3r 2 )mi A /2c 2 

q 2 “ (7-2PA/c + r 1 r 2 i 2 /4cm)£ 2 m/c 

q 3 “ (r j + 10r 2 ~ 3r 2 PX/c)£ 2 /2c 

q 4 = 1 

Routh-Hurwitz criteria lead to the critical load 
^md " 35 / 12 + M*^/6p<2 + + 

^J^5/12 + ^i^/6^ 2 - p» 2 (l + + 3p^p 2 )/4p,J 

+ (^ x + 3p, 2 ) 2 /6n 1 n 2 j- 

where 


(4.144) 


(4.145) 


(4.146) 


(4.147) 


(4.148) 
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F - PjJ/c 


(4.149) 


jij - rjje/(4cm) 1 ^ 2 (j ** 1,2) 

A point of interest regarding (4*148) is that although F^ -» » as -» «> 

(i.e., the system can be made stable for arbitrarily large P), it remains 
finite for jj,^ / 0 and ^ In fact 


lim F 
^2 


md 

oo 


f 10/3 if ^ + 0 is finite 
(10 + X)/3 if p,j ** 


(4.150) 


becomes 


In the case of small damping, i.e., p,,. « 1, F^ is written F^ and 


md 


F md ° 35/12 + h/6 + 1/4h 

- [35/12 + h/6 + 1/4h) 2 - (4k + 1/k + 45)/6] 1/2 (4.151) 

where 

K - Hj/^2 

We note that the greatest destabilizing effect is realized as h -> 0. We further 
note that F^ is a monotone increasing function of h. 

Comparison with the case of internal viscous damping (Eq. 4.14) reveals 
that becomes unbounded as either or becomes large, provided the other 

parameter is non-zero. With magnetic damping, on the other hand, we have the 
result (4.150), and therefore magnetic damping can be said to be weaker than 
internal viscous damping. Furthermore, while it is possible to eliminate the 
destabilizing effect with a viscous damping coefficient ratio of 11.07, the 
critical load of the magnetically damped system is always smaller than F^. 

It may be of interest to compare the effects of magnetic damping with 
those of linear external viscous damping. If in the double pendulum system 
of Fig. 4.24 external damping forces act which are proportional to the velocity 
with constants (along OA) and k^ (along AB), then the damping force distri- 
butions are, for small angles, linear functions of distance along the rods (see 
Fig. 4.25). In this case, the equations of motion are 

3 m£ 2 cp^ + (kj+ 3k£) q>j£ 3 / 3 - (P£-2c)cp^ 

+ me 2 jL + kiJ 3 /2 + (P-e-c)cp- = 0 (4.153) 

cont. 
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(4.153) 


ml? cpj + il^/2 - ccpj^ + raJL^ $2 + ^2^2 J ^ 3 + Cc P2 ” 0 

A development parallel to that which led to Zq. (4.148) yields the fol- 
lowing expression for the critical load parameter: 


F ey - (L - J L 2 + 4HN ) /N 


(4.154) 


where 

L - 4(T1 X + 8T1 2 )(T1 1 - 21^) + 5T) 2 (T1 1 + 3T1 2 ) (7 + 4^71.,+ 3 ^) 
M - 2(T| 1 + 8T| 2 ) 2 + Oljd- 3T) 2 ) 2 

- (\+ 8 V (\ + 3 V (7 + 4 %V 3 \ 2) 

N - 10T[ 2 (2T1 1 + T) 2 ) 

and 


(4.155) 


Tlj » k^ 2 /(36cm) 1/2 (4.156) 

Examination of (4*154) establishes the following results: 

Just as in the two previous cases, the system can be made stable for 
arbitrarily large P by letting be arbitrarily large, i.e* f as -» » 

(0 j* 7] 2 Is finite) -> The behavior of (4.154) as both T)^ and T| 2 become 

large resembles that of rather than F^ v * In this case 


lim F 

ev 

\ 00 


r 16/5 if T1 1 
(16 + 2\)/5 


is finite 

If \ - X7} 2 


(4.157) 


Thus the external viscous damping of the type being considered is also weaker 
than the internal damping. A It is noted, however, that when both damping pa- 
rameters become unbounded F - < F for all X >4. 

ma ev 

Another feature common to all three types of damping is that the critical 
load approaches the value 2 as ^(I^,^) approaches zero, and this result is 

independent of ^(Bp Tlj) • If, however, |jl^ "^1^ a PP roaclies zaro * both F^ and 

approach the value 1/3 (independently of pt 2 > ® 2 ) while tlie value of 

depends upon T| 2 r 
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(4.158) 


and 


lim I = 3(7 + 3TL 2 )/2 - 32/5 - { [3(7 + 3T) 2 )/2 - 32/5 
*n 1 ~» o ev i. l ^ 

+ 274/5 - 48(7 + 3T| 2 2 ) /sj- 1/2 


lim 

V° 


/lim F 


ev 


) 


2 


lim 

ti 2 -» ® 


/lim F 
0 


ev 


) = 16/5 


(4.159) 


In the case of small damping, i.e.. 


Tlj « 1, 


(4.154) becomes 


F ev “ C35 (k + 3) + 4(h + 8)(k - 2)]/10(2k + 1) 

- { + 3)/2(2k + 1) + 2(k + 8)(h - 2)/5(2h + l)J^ 

r 9 ”1 1 1/2 

- ^8h + 78k + 62 J /5(2* + 1) j (4.160) 

where 

H = \/\ (4.161) 


This result differs markedly from the cases of internal viscous and magnetic 
damping in that F does not depend upon h but is equal to the constant value 

A 

of 2. Therefore, the critical load parameter F^ in the case of external vis- 
cous damping differs from the value of 2 at most by terms which are of second 
degree in 7^ . 


4.6.3 Magnetic Damping in a Continuous System 


As a second example of the effect of magnetic damping, an elastic conti- 
nuous cantilever acted upon by a follower force P will be considered (cf. Ref. 
[41]). According to the Bernoulli-Euler theory, the equation of motion which 
describes the system when El is constant is 


*4 *2 .2 

WT O V . 0 V . o V , X 

El 7 + P — o + p — i “ v(x) 




3x‘ 


at* 


(4.162) 


where El is the flexural rigidity, p the lineal mass density, and w(x) is the 
force per unit length acting in the y-direction. 
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The displacement v(x, t) must satisfy the boundary conditions 
V(°,t) “ -g (0,t) - 0 

<i,t) U,t) ■> o 


(4.163) 


Now, if the cantilever is a portion of an electrical circuit having re- 
sistance R and if the system undergoes motion in the presence of a uniform 
magnetic field $ o whose direction is normal to the x-y plane, then the in- 
duced current will provide the following damping force distribution: 


w(x) 




dx 


(4.164) 


With the substitution of (4.164) into (4.162) and the introduction of the 
dimensionless parameters 

g - x/j£ 

T - t(pA 4 /EI )" 1/2 

(4.165) 

F - Pi 2 /El 

V = $ o V/(pR 2 EI ) 1/2 


the equation of motion and boundary conditions appear as 

,1 


^ + F^ + „ 

ar ar 


J «« 


d 2 V 

St 2 


v 



0 at | - 0 


(4.166) 


a 2 v 


a 3 v 


ar 


0 at g 


In order to deduce stability criteria for this system, ve consider modal 
solutions of (4.166), i.e., we set 


v(g, t) - Y(g)e' 


Ot 


(4.167) 
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Substitution of (4.167) into (4.166) results in the following boundary value 
problem for Y. 


d 4 Y , - d 2 Y 
4 * 2 


Y 





- Q 2 (K-Y) » 0 
at 5*0 
0 at ? - 1 


where the functional K has been defined according to 


Ov 



o 2 k 


(4.168) 


(4.169) 


(4.170) 


We proceed now in a purely formal manner to solve (4.168) subject to 
(4.169): The general solution of (4.168) is 

Y ° A^ sin + A^ cos + A^ sinh + \ cosh u^g + K (4.171) 


where 

m 2 - [(F 2 - 4Q 2 ) 1/2 + F]/2 
(0 2 2 - [(F 2 - dn 2 ) 1/2 - F]/2 


(4.172) 


Substitution of the solution (4*171) into the boundary conditions leads to a 
system of four nonhomogeneous algebraic equations in the coefficients 
A. ( j “ 1,*..,4) whose solution is found to be 

2 3 

Aj ■ - KuUj u>2 (tu-^sin ujjcosh <u 2 + <« 2 cos tu^sinh m 2 )/A 

3 2 2 

A, » Kuj. to, (ou.ou-sin (D.sinh co,- U). cos ui. cosh uu- od, )/A 
i l z l z i i z (4.173) 

3 2 

Aj ■ Kuj^ u> 2 (cousin cohosh u> 2 + (« 2 cos (DjSinh <» 2 )/a 

3 2 2 

A^ *■ - Ka>^ cu 2 (cousin to^sinh ou 2 + <« 2 cos (UjCosh u) 2 + tu^ )/A 

where 
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A ■ + 0)2 + 20 ^ cu^ cos ou^ cosh 

+ ^ 2 2)8in sinh (4.174) 

The fact that the above solution is given in terms of the unknown func- 
tional K(Y) is not a severe deficiency since primary interest in the present 
context is focused on the nature of the complex frequencies ft. As long as 
the real part of ft is negative, the rod will be asymptotically stable, i.e., 
will oscillate with exponentially decreasing amplitude. The characteristic 
equation, by which the nature of ft may be examined, is obtained by requiring 
that (4.171) satisfy (4.170) nontrivially, the Aj being given by (4.173). 

This requirement leads to the following transcendental equation for ft: 

T 4 2 3 

Aft * v I i (1-cos u)j)8inh sin cosh 

3 2 4 

- («2 cos ou^ sinh sin (^(l-coshu^) 


- 


sinh 




- 0 ) o sin UD 


j 


vA„ 


(4.175) 


For small values of F and for positive damping (v > 0), all the roots of 
(4.175) are located in the left half of the complex plane. As F is increased, 
one of the roots approaches the imaginary axis and subsequently takes on a 
positive real part. When this occurs then, by virtue of (4.167), oscillations 
with exponentially increasing amplitude will result. The value of F beyond 
which ft has a positive real part will be designated as F^. As F -» F^, there 

must be one imaginary root of (4.175), and since are real whenever ft 

is imaginary, both sides of (4.175) must approach zero simultaneously. 

The critical load is found numerically to be F ^ * 12.84, and the result 

is independent of the magnitude of the magnetic damping (provided it is non- 
zero) • Comparing this value with the critical load in the absence of a mag- 
netic field, F = 20.05, we find that the magnetic damping has a destabilizing 
ex 

effect of 36 percent. 

4.6.4 Retarded Follower Force 


The system with two degrees of freedom discussed in Sect. 3.1 was sub- 
jected in Ref. [42] to a retarded follower force. It is remarkable that a de- 
stabilizing effect is associated also with retardation of a follower force with 
constant time lag t, which was specified as 

0(t) - cp 2 (t-T) (4.176) 
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The linearized equations for small motions about the position of static equi- 
librium (cp^ *= cp£ * 0) are 


,2 - 


3m i cp^ + ml + (2c-P/)cp^ 


cc^ + P£0 * 0 


(4.177) 


ml* cp^ + ml* cp^ - ccp^ + (c-Pi)cp2 + P IQ *= 0 

Solutions are sought again in the form (4.4) and lead to the frequency equa- 
tion 


3ml* w* + 2c - VI 

,2 2 

mi CD - c 


m£ 2 u> 2 - c + PXe" WT 
mS, 2 w 2 + c - P i + Pj ee _a,T 


- 0 


(4.178) 


The presence of exponential terms suggests the application of Pontryagin's 
stability criteria rather than those of Routh-Hurwitz. After a comprehensive 
and rather elaborate analysis we arrive to the important result that a very 
small (vanishing) time lag renders the system unstable for all positive (com- 
pressive) values of the applied force P. Even under the most favorable time 
lag the critical load was found to be F = Pi/c ** 0.177, as compared to 
F * 2.086 for the same system without any retardation t. 

In supplementing the analysis of Ref. [42], it may be remarked here that 
for small time lag the stability investigation can be readily carried out em- 
ploying the simpler Routh-Hurwitz criteria. If ^(t-r) is expanded into a 
Taylor series about cp^(t) (see Ref. [18]) and if only the first two terms are 
retained, the equations of motion simplify to: 


3m£ 2 + ml* cp 2 4* (2c-P£)cp 1 - ccp^ + P£<P 2 ~ “ 0 

(4.179) 

2 „ 2 ~ 

mi, cp^ + ;p 2 - ccp^ + (c-PjOep 2 + = 0 


and lead to the frequency equation of the form 
P 0 Q 4 + PjO 3 + p 2 q 2 + p 3 n + p 4 = 0 

with 

Pq = 2; Pi » - 2FT; 

P 2 = 7 - 2F; P 3 = FT(F-3); p 4 * 1 


(4.180) 


(4.181) 
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where 


Cl 2 *= ml 2 wta F - P£/c; 
T = TU>/ft = t /c/£ /m 


(4.182) 


One of the Routh-Hurwitz conditions for asymptotic stability is that p Q and 

p^ be of the same sign which results in the critical value of F cj _ * 0, that 

is for F < 0 the system is asymptotically stable and for F > 0 it is unstable, 
verifying thus the result of Ref. [42]. The remaining Routh-Hurwitz condi- 
tions do not supply more stringent requirements of F. It is noteworthy that 
F^^ does not depend on the value of T. The conclusion Is reached that in the 

presence of even the slightest lag the system is unstable under a compressive 
follower force. Further, the Taylor series expansion introduced above clearly 
exhibits that small time lag is associated with the introduction of terms of 
odd power in the frequency equation, having this in common with linear viscous 
damping. 

The destabilizing effect is in general introduced by any sufficiently 
small, velocity-dependent forces, such as, for example, Coriolis forces. Some 
types of such forces are realized by fluid jets; they have received conside- 
rable attention and have been referred to as "jet damping," [43]. Others are 
produced, e.g., by flow through pipes [44]. 


4.7 Uncertainties 


The foregoing examples of various destabilizing effects amply illustrate 
the necessity of a firmer grasp of certain aspects in the analysis of stability 
problems as applied to systems subjected to follower forces, cf. [45]. What 
is needed in particular is additional insight into the experimental determi- 
nation of system parameters, cf. Sect. 8.2. If very small, even vanishing 
quantities which induce a destabilizing effect have such a decisive influence 
on the critical loads calculated analytically, how should these quantities be 
measured with required accuracy? Further, how can one be sure that the "correct" 
or "right"parameters have been included? It is even conceivable that (vanishing) 
destabilizing effects exist which have never been thought of as yet, and which 
have perhaps an even stronger influence on the stability boundaries of a given 
system than any of those mentioned. Below some attempts are described to re- 
medy this obviously unsatisfactory state of affairs. 

If we are dealing with a man-made system, it would probably be desirable, 
if possible, to make it well-behaved by means of a suitable choice of system 
parameters and, in particular, by making it strongly asymptotically stable to 
begin with. This is done sometimes in control systems where the "doubtful," 
"critical" or "marginal" case of Liapunov (pure imaginary roots of the charac- 
teristic equation) is interpreted as describing an inherently unstable system. 

It is well known that if a system is asymptotically stable, small "destabi- 
lizing" quantities will have but a small effect on the critical loads; this 
effect will vanish with the vanishing of the "destabilizing" quantity. 


73 



In many man-made systems and in given natural systems the uncertainty can- 
not be circumvented in this manner. It may then be suggested that the analysis 
of stability be replaced (or supplemented) by an analysis of "patterns of be- 
havior" of the disturbed system for various ranges of the controlling parameter 
(force). Since we are interested here only in oscillatory response to a dis- 
turbance, three types of behavior are qualitatively sketched in Fig. 4.26. Let 
it be our aim to classify the response in just two categories. Depending upon 
the specific performance requirements of the system at hand, it may be meaning- 
ful to place the response in Figs. 4.26a and 4.26b into one category and the 
response in Fig. 4.26c into the other. In the first category the disturbance 
remains small during a certain interval of time, while growing fairly large in 
the second category during the same interval. 

As a measure of the rate of growth of the oscillations it is convenient to 
introduce the largest real part a of the relevant root of the characteristic 
equation. This is analogous to the introduction of the smallest negative part 
as the "absolute stability margin" [18]. 

The stability analysis corresponding to a given or > 0 can be carried out 
by introducing the transformation (Fig. 4.27) 

0 ■ p + a (4.183) 

into the characteristic equation, e.g. , 

P o n 4 + Pjft 3 + P 2 ft 2 + P3^ + P4 53 0 

which yields the modified characteristic equation for p 

a Q p 4 + a^p 3 + a^p 2 + a 3 p + a^ = 0 (4.184) 

where 

2 

a o = p o ; a i 38 p i + 4<*P o ; a 2 = P 2 + ^l + 6of P o (4.185) 

a 3 “ P 3 + 2orp 2 + 3of 2 Pj+ 4or 3p o ; a^ = p 4 + orp 3 + <* 2 p 2 + 0 ?^+ <x\ Q 

Applying the usual Routh-Hurwitz criteria to this modified equation, the cri- 
tical force can be calculated. For the system of Fig. 4.1 the critical force 

F is to be calculated from 
c 

2 2 

x ■ a^a^a^ - a^a^ - a^ a^ (4.186) 

and is found to be 

F c = (B ( t J (4.187) 
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where 


and 


A * 8or(d 1 -2ad o ) 


B - d x d 3 + 2«d 1 d 2 - 4«d o d 3 - a 2 d* 

c - W3 - d 0 d 3 2 - d 3 2 d 4 


(4.188) 


2 

d o m P o ; d i * Pi + 4ap o ; d 2 = 7 + 8^2+ 30^+ 6a P Q 

d 3 - p 3 + 2a(7+B 1 B 2 ) + 3a 2 ?! + 4a 3 p Q (4.189) 

d 4 “ P 4 + «P 3 + or 2 (7+B 1 B 2 ) + + a 4 p Q 


The results of the numerical calculations are displayed in Fig. 4.28. Thin 
solid lines represent the critical force F c as function of the growth para- 
meter or for given damping coefficients B^. The thin curves in Fig. 4.28 are 

the same as those in Fig. 4.5, but they have been calculated in a different 
manner and their interpretation is also entirely different. 

The critical force as defined with the aid of the growth parameter a is 
not entirely satisfactory because it does not separate the different types of 
behavior illustrated in Fig. 4.26. We seek now to define what we may call a 
"transition” force F fc below which a (whether positive or negative) would be 

relatively small and above which it would be relatively large. It appears to 
be reasonable to define the transition force F fc as the force for which the 

absolute value of the curvature of a given curve F(<*) = 0 attains a maximum. 
Corresponding calculations have been carried out and the values of the tran- 
sition force F fc for various values of the damping coefficients have been 

joined by a thick solid line in Fig. 4.28. It Is noted that for given damping 
coefficients there exists an associated transition force which in turn 

corresponds to a certain particular value of growth parameter Qf t . The expe- 
rimental determination of system parameters associated with F fc appears to be 
feasible. It should be also observed that as damping decreases F £ approaches 
F fif while approaches F^ as damping increases. The relationship between 
critical and transition forces is thus clarified. 


It is rather evident that the two types of system behavior which are 
separated by F fc can be less or more different and thus It may be appropriate 

to introduce the notion of degree of separation ct associated with any particu- 
lar value of F t * This additional characteristic separating "patterns of 
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behavior" could be made to depend on the magnitudes of curvature and slope of 
the function F(or) at a^. 

The parameter a can be employed yet for another, purely mathematical 
purpose. For small of, and Eq. (4.186) may be written in the form 

[or 2 + a(B 1 +6B 2 )/4](F-F el )(F-F e2 ) - (B^) (B^B^ (F-F d ) =0 (4.190) 

Here F^ and F^ are the critical loads obtained in the absence of damping, 

while F^ is the critical load for vanishing damping. It is observed that either 

critical load may be obtained by a limiting process in Eq. (4.190), which 
gives either o', F, or B^ in terms of the remaining three quantities. If 

the growth parameter a — 0, it is seen that F F d , which depends only on 
B 3 B ]/ B 2 but n0t ° n B 1 an< * B 2 contrast, if the damping coeffi- 

cients are made to vanish first, then F-*F- orF-*F 0 , regardless of the 

ei e2 

value of the (small) value of or. Thus the introduction of the growth para- 
meter a permits to approach the critical load for no damping F^ even in the 

presence of vanishing damping, removing mathematically any destabilizing 
effects. 
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CHAPTER V 


CONTINUOUS SYSTEMS 


5.1 Introduction 

In the preceding chapter one aspect of problems of continua was not ela- 
borated upon, namely, that stability must necessarily be defined with respect 
to a metric (sometimes implied) which measures distance in an infinite- 
dimensional space. One has to define what is meant by ,f nearness n to the 
equilibrium configuration whose stability is being examined. This metric 
may be postulated in various suitable forms, depending upon the physical as- 
pects and the requirements of the specific problem at hand. The equations of 
the boundary value problem of a continuum, together with an explicitly defined 
metric p form a functional metric space whose fundamental properties depend 
strongly on p and thus lead to different results of a stability analysis. 

With reference to a conservative system, Koiter [46,47] has pointed out 
that a conventional generalization of Liapunov* s definition of stability, 
which requires that the displacements and the velocities remain arbitrarily 
small at each point and for all positive time, provided the initial disturb- 
ances are sufficiently small, can hardly be considered satisfactory. A modi- 
fied concept of stability was suggested which reduces to Liapunov's definition 
for the case of a discrete system. 

In this Chapter, following the development of Ref. [48], a sufficient 
condition for the stability of a linearly viscoelastic continuum subjected to 
surface tractions which follow partially the deformation of the solid is estab- 
lished with respect to an average metric. 


5.2 Definitions of Stability 

We consider a finite isotropic, homogeneous, linearly viscoelastic solid, 
bounded by a regular surface S, contained in a volume V. At the time t = 0, 
the solid is in a state of initial stress ct^: i,j * 1,2,3, caused by a sys- 

tem of partial follower surface tractions p^, applied at the boundary S. We 

shall refer to the state of initial stress of the solid as unperturbed (equi- 
librium) state and study its possible motions with reference to this state. 
Furthermore, we shall assume that the quantities describing the perturbed state 
are small (these quantities will, subsequently, be indicated by a bar) so that 
all terms of order higher than the second may be neglected. The equations of 
motion of the perturbed solid, referred to a fixed orthogonal Cartesian coordi- 
nate system, are [8] 

\j,j + (CT jk^i,k>,j - m “i " 0 in v > ^ij n j + a jk“i,k n j “ ?! on s > * 

i,j,k « 1,2,3 
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where m is the mass density, x^ are the coordinates, the displacement com- 
ponents measured from the unperturbed state, n^ the components of the unit 
normal to S, p^ the perturbations of the applied surface tractions. In these 

equations and in the sequel the repeated indices are summed over the range of 
their definition. A comma followed by indices k, j indicates differentiation 
with respect to x^, x^, and dots denote derivatives with respect to time. We 

shall assume here that 




on 


(5.2) 


where a(x) = cKx^x^x^) is a parameter which serves to describe the manner 

in which the surface tractions follow the deformation. If a = 0 the system 
is conservative and for a = 1 we have the case of follower force introduced 
in [8]. The constitutive equations shall be taken in the form 

a ij * C ijki u (£,k) + C (ijk£) u (£,k) * °(k,£) ~ 2 


'ijk* - X6 1:j 6 w + 2*6. k 6 jje 


(5.3) 


C.' . = X'6. .6, , + 2p, / 6., 6.„ 

ijkjfc xj k l p lk jX 

where 6^ is the Kronecker delta, \ and p, are Lame constants, and \ i and p/ 
are viscous constants corresponding to Lame constants. 

A general solution to the nonself-adjoint mixed initial and boundary 
value problem (5.1) cannot, in general, be easily obtained. Therefore, in 
order to study the stability of this system, we have to resort to some other 
means and, consequently, we shall not expect to gain as much information con- 
cerning stability as we would if we were to construct and evaluate a general 
solution of the system. As we shall see, this is by no means a shortcoming. 

A strong stability criterion,, that may be imposed on the system and which 
could be applied if we were to solve system (5.1) completely, would be of 
doubtful interest. 


In this connection, we shall consider a certain functional (which, in 
effect, expresses the energy of the system) and explore the stability of 
(5.1) in some appropriate average sense. Furthermore, we shall show that 
the usual Galerkin method, which reduces the system of partial differential 
equations (5.1) to a set of ordinary differential equations, yields the same 
results as those obtained by a study of the functional mentioned, provided 
all the series expansions employed converge in an average sense. 


To this end, we consider a complete set of normalized eigenvectors, ob- 
tained by solving the homogeneous, self-adjoint system deduced from (5.1) by 
~ =r' = « « o, which has the same geometrical boundary 


setting a 


ij ijkX 


= P 
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conditions as the original problem. Let this set of orthonormal eigenvectors 
be denoted by {cp in (x)]; i “ 1,2,3, n * 1,2,...,®. We shall reduce our origi- 
nal system of partial to a system of ordinary differential equations by ex- 
panding Uj, and its derivatives in terms of these eigenvectors, without any 

attempt to resolve the question of convergence. However, some comparison 
between the results obtained by applying this method to some simple problems 
and the exact solutions [8] certainly suggests that convergence may be as- 
sumed.* In our problem, we shall therefore state that if convergence exists 
(in an average sense at least) then the two methods yield Identical results. 

Let us now consider the fundamental question concerning stability of a 
solid, and review first the definitions of stability for a discrete system, 
mentioned in the Introduction. 


We examine a system with r degrees of freedom described by generalized 
coordinates q^ and generalized velocities q^; n * 1,2, ...,r. For a holonomic 

and autonomous system, we write the equations of motion as 

z n “ f l* ^2 * * * * ^2r^ * ^ B l,2,...,2r (5.4) 

where 


z 

n 


q 


n 


z 


r+n 


v 


n 


1.2 


, • • • , 


r 


and f R (z) are bounded, continuous, real functions vanishing for * 0, We 

assume f n satisfy all the conditions required for the existence of a single- 
valued solution for t > 0 in the region of the definition of z r . Furthermore, 

we represent the state of this dynamic system by a point in a 2r-dimensional 
Euclidean space, E 2r* wit ^ coordinates z^; n *= l,2,...,2r. The equilibrium 

state of the system at the origin is said to be stable if for any c > 0 we 


2r 


can find a 6 > 0 depending on c only such that when ) z ^ < 6 at t c 0, 

_ L-i n 

2r n=l 


V 2 

we have ^ z n < * for all t > 0. In the opposite case z q = 0 is called 
n=l 

unstable [14]. Furthermore, z * 0 is called asymptotically stable if it is 
2r n 


stable and lim 

t -» ® 


K 2 


n*=l 


0 . 


* The paradox in the problem of flutter of a membrane, as was shown in 
[8], is not related to the fact that the system is nonself -adjoint. 
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The above definitions of stability are due to Liapunov [11,14]. He also 
supplied the proofs of necessity and sufficiency, employing the notion of dis- 
tance in the finite -dimensional Euclidean space E 2 r » 

For systems with an infinite number of degrees of freedom (continuous 
systems) the notion of distance in an infinite dimensional space needs to be 
introduced, if one wishes to extend Liapunov's concepts to such systems. In 
this case, we have to be concerned with functionals rather than functions and 
must explicitly define a measure (metric p) of distance of two states of the 
system and then study the stability of the system with respect to this metric p. 
The metric p may be selected in any suitable manner (provided it satisfies three 
fundamental conditions [49] so as to fulfill some physical requirements of the 
problem at hand. It may be desirable, for example, to limit the displacements 
and the velocities at each point of the solid, in which case we define 

. . 

= ^ i u i + u i u i everywhere in V and on S 

In some other cases, we may wish to restrict the strains as well as the dis- 
placements and the velocities at each point of the solid, such that 

P 2 = + u i u i + u i j u i j everywhere in V and on S 

For most practical problems, however, it is usually preferable to define 
p in an average sense; for example, 

p 3 = / v + + z ±\] dv 

We now state the definition of the stability of the initial state of a 
solid with respect to an explicitly defined metric p, by appropriately ex- 
tending the corresponding definition for a finite system. 

The initial state of the continuous solid is said to be stable if for a 
given e > 0 we can find a 6 > 0 depending on e only such that when p < 6 at 
t = 0 we have p < e for all t > 0. In the opposite case, the initial state 
is called unstable. Furthermore, the unperturbed state is called asymptoti- 
cally stable if it is stable and lim p = 0. The sufficiency theorem of stabi- 

t 00 

lity may now be stated as follows: 

Theorem , In order that the unperturbed state of system (5.1) be stable 
with respect to a metric p, it is sufficient that there exists, by virtue of 
the requirements of the boundary value problem (5.1), a finite, nonincreasing 
functional which is identically equal to zero for p = 0 and admits an infini- 
tely small tipper bound with respect to the metric p. 

This theorem is an appropriate version of the theorem of stability given 
by Movchan [50]. In the sequel we shall use this theorem to establish a suf- 
ficiency criterion for the stability of system (5.1). But let us first dis- 
cuss some aspects of the definition of stability. 
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It is seen that the stability criteria are highly dependent upon the spe- 
cification of the metric p. We may not, therefore, expect to apply a criterion 

obtained, say, for to p^ and get like results. The problem which was treated 

by Shield and Green [51] may exemplify this very point. An isotropic, homoge- 
neous, linearly elastic sphere was perturbed by radially symmetric applied in- 
finitesimal disturbances at t B 0 and it was shown that the strain at the center 

of the sphere can become finite for some t > 0. Let us show that although this 

system is unstable with respect to the metric p^> it is stable with respect to 

p^. To this end consider the following functional 

H 1 " 2 [ J v <“i“i + dv ] 

whose time derivative is zero by virtue of the equation of motion, and which 
admits an infinitesimal upper bound with respect to the metric p~. From the 
inequalities [52] 

c i J T Vi dv s J„ “i.j”i.J dv 
C 2 J y “i,j”l,j dv !„ °lJM“l,A.* dv 

which are valid for all admissible motions of the solid, with and being 

fixed positive constants independent of u^, we immediately construct the in- 
equality 


£ Kp^ for all t SO 

where K is also a fixed positive number not dependent on u^. We let < Ke 
and obtain p^ < e at t 3 0. But is a nonincreasing function of time. 
Therefore Ke is an upper bound of for all t £ 0, which implies 

p^ < e for all t ^ 0 

In [51], the initial disturbances were taken to be 

u “r“fr = 0, “ = r^ [r at t = 0 

where r measures distances from the center of the sphere, c * — , and 

f(r) is given by 

f(r) ■ 0 0 £ r £ a 

1 4 

f( r ) = ^ ■ (r-a) (r-a-2ea) a £ r £ a + 2ea 

e a 

f(r) =0 a + 2ea £ r 
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A simple calculation shows that * 0(e) at t » 0, Furthermore, at t ■ a/c 
we have, for 0 £ r £ 2«a, 

u - -y-j r 2 (2ea-r) 3 (7r-6«a) 

« a 


which immediately yields “ 0(c) at t * a/c, while the strain at the center 
of the sphere at this instant is finite: 



r*ea 


1 , 



6 


In this example, one is able to obtain an exact solution to the differen- 
tial equations of the boundary value problem. Therefore, one is in the position 
of requiring as strong a stability criterion as one pleases. We see that the 
system is not stable with respect to although it is stable with respect 

to py The important point to note in this connection is that the stability 

with respect to the metric p^ could have been deduced without possessing an 

explicit solution of the problem. 

In most practical problems, the system may well be stable for all practi- 
cal purposes, while it may not satisfy the pointwise stability conditions with 
respect to the metrics p^ and p^. In those cases there may exist a finite 

number of points in V where an infinitesimal perturbation at t ■ 0 may cause 
finite, say, strains at these points for some t > 0. If the collection of 
these points forms a set with measure zero, then the stability may exist with 
respect to the metric p^. 

The metric p^ seems to be more appealing also from a purely mathematical 

point of view. In this regard, let us note that the series expansion of a 
piecewise continuous function in a finite domain is an approximation in a mean 
square sense and not a pointwise representation. The following discussion will, 
therefore, be devoted to the stability of system (5.1) with respect to the 
metric p^. 


5.3 Analysis of Stability 

We consider a functional H given by 

H ■ H I s [ "Vi + 

+ 2 l j"k,i ' 

and note that, from the requirements of the 


+ <l^)o Jk « ltk u lfJ ] dv 
°jk WB i.j ) .k«i3 dv dt } 

boundary value problem 


(5.5) 

(5.1), H is 
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a continuous functional which vanishes identically at the Initial unperturbed 
state of the solid, » 0. The total time derivative of H is 


dH 

dt 


J>Vi + c ijwA,A,t + w Vi,A,] 


But we have 


+ C ijkA,j”k,i - V’A.AkV dv 
J 0 [c iik«A,k + 'iw'vi 1 A,j iv - L va® • J„ 5 u,j“i dv 


(5.6) 


and 


J v { <l - a>a jk“i,A,j - a jk<« u i,j> f k ; i} dv - ■ J v ta jk“i,k^A dv 

+ J s [ff jk“i,k n j ■ ^i.k^ V s 

where in the last reduction we have used the fact that for the unperturbed 
state we have 

c ij, j “ 0 in V and CT jk n k “ p j on S 

Equation (5*5) now becomes 

S? " I v {A - j - (CT jk“i,k>,j} V v + J s t5 ij n j + a jk“i,k n j 

- u^S (5.7) 

which is identically equal to zero by virtue of Eqs. (5.1) for all actual per- 
turbed motions of the solid. Moreover, if H is a positive definite functional, 
then it admits an infinitely small upper bound with respect to p^. To show 

this we let |u i | < /e, | | < /* and ju^ ^| < /c at t » 0*, and obtain 

p^ < 15Vc at t * 0 

Then, as H > 0, we have 

HCKe-6 at t * 0 


* The initial disturbances may also admit singularities 
of isolated points in V such that » 0(c), and H *= 0(e) at 


at finite number 
t » 0. 
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where K is a positive constant. But 6 is an upper bound of H for all t > 0, 
as H is a nonincreasing function of time. Therefore, if H is a positive de- 
finite functional, then all the requirements of the sufficiency theorem are 
fulfilled and we have the following theorem: 

Theorem . For a linearly viscoelastic solid subjected to a set of par- 
tial follower forces to be stable with respect to the metric p^, it is suf- 
ficient that the functional H given by Eq. (5.5) be a positive <Ief ini te quan- 
tity for admissible perturbed motions of the solid about the state of initial 
stress. 

Let us note that the requirement of H being a positive definite functi- 
onal may imply a stronger stability condition than is given by p^. This 

touches then upon the question of the necessary conditions which will not be 
dealt with here. 

From the above discussion we may conclude that the commonly used energy 
methods yield stability criteria with respect to an average metric p^. There- 
fore we may not, by any means, expect to retrieve any more information than is 
retained after this averaging process. This conclusion is also valid for most 
approximate methods such as the Ritz, the Galerkin, and other methods, where 
we use some averaging processes to reduce the system of partial to a set of 
ordinary differential equations. We shall explore this point further in the 
sequel, but let us make first another remark regarding system (5.1) and func- 
tional H. We let solution of (5.1) be of a form u^ - ^(x)e pt and obtain from 
(5.5) 

K ' ,2 ' >t { 5 + + < 1 -“>Vi.k*i. J 3 dv 

+ dv } (5 ’ 8: 


If we substitute u^ = \|f^e pt into Eqs. (5.1), we obtain an eigenvalue problem 

with eigenvalues p. From Eq. (5.8) we may conclude that, for H to be a non- 
increasing function of time, p must have a nonpositive real part. 


We now reduce Eqs. (5.1) to a set of ordinary differential equations. We 
assume that u and its derivatives can be expanded in terms of the complete 

set of eigenvectors fcP in ( x )} * “ 1,2,3, n « 1,2,...,«, such that 


i 


N 

I Vi' I c Pin<Pin q n 2 < t) 

n=l 


dv < e., 


L U i U i'I c Pin C Pin q n 2(t) 


dv < «„ 


n=l 


(5.9) 

cont. 


84 


N N 


*V 

and 


dv < Co 


J | U i t j U k,i" Z Z ^in, j Cp km,i q n^ t ^ q m^ t ^ 

V n=l m=l 

N N 

1 K.jA ' ^ ^ C Pin, j k c Pira q n (t >%< t) | dv < *, 


(5.9) 


n *=l m»l 


N N 


J v K.A.i- Z I 'Pin,j'fW,X q n( t > q m( t >| dv < «5 * 1 ’ 2 ’ 3 

n=l m=l 

for some N > M, where M is a large positive number depending on e i ; i*l,2,...,5 
in the above inequalities and may be made as small as we please by selec- 
ting M sufficiently large. For such an M, Eq. (5,7) reduces to 


N 

i 

m=l 


I {%n + l °mn q n + % 1 ( V b mn> q n } Sn = 0 


n=l 


n=l 


where 


b mn = 72 [ I v ( 1 - a)a jk c Pin > j t Pim,k dv ' J/jk^in. j> f k*im dv ] (5 * 10) 
m 


and 


C mn J v C ijkX qJ kn f / q> im, j dv 

In obtaining (5.10 ) 9 in addition to the Gauss theorem we have also utilized 
the fact that {cp in } are solutions to 

C ijk*'Pkn„ej + mU, n 2c Pin “ 0 in V * C ijki^n,je n j = 0 on S ’ J y m 'Pin c Pim dv " 6 , 

For ^ ; m = 1,2,...,N not identically zero, Eqs. (5.10) yield 


N 


5. + 1 + >% 2 1 <W » 

n=l n“l 


(5.11) 


which is a system of nonself-ad joint, ordinary differential equations. 
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Similarly, H reduces to 


N 

i 

m=l 


N 


“"I I { [C + ».V] + I UV 2 J 1 [ G mAflm + l wM^\ dt } 


n«l 


H 

z 

n*l 


where 


a mn “ I ( 1 - Qf)CT jk < Pin,j c Pim,k dv * b mn * " J/jk^in, 


and 


b 

mn 



(a + b ) 
mn mn 


(5.12) 


For a positive definite H in a region < R; R > 0, we can find an M 

such that H is also a positive definite quantity within a ring < *p^ < R > 
where “£3 is defined by 

*3 - Z <’n 2 + C> 

n=l 


in a 2 N-dimensional Euclidean space. Moreover, R^ is dependent only upon 

in inequalities (5.9) and may be made as small as we please by choosing M 

large enough. From the stability theorem we therefore conclude that, for 
system (5.1) to be stable with respect to the metric p^, it is sufficieit 

that H be a positive definite quantity. But H vanishes for "p^ * 0 and dH/dt 

is identically equal to zero along any path satisfying equations (5.11). 
Therefore, by Liapunov's stability theorem [14], system (5.11) is stable 
when H is a positive definite quantity, and likewise when Q is a positive 
definite quantity. 


The study of stability of the system of linear homogeneous ordinary dif- 
ferential equations (5.10) is, however, a classical mathematical problem. 

For the stability of (5.11), it is necessary and sufficient that the roots 
of the characteristic equation of _(5.11) have nonpositive real parts. How- 
ever, the study of the functional H, which in fact is a statement of the 
energy of the system, can provide us with a better insight into the physical 
behavior of the system. Divergent motion may occur if, for a virtual (static) 
displacement of the system, the work of the applied forces equals the change 
in the strain energy of the system, namely. 


4 I 2 ['ukt-M-M + dv ‘ I < ’jk ( “ 


“mVv” 
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or equivalently 


6 J v 2 tc ijkje“i,j“k,je + CT jk“i,k“i,^ dv - J s “Pj“i,j 6 V s “ 0 < 5 - 13 > 

where 6 is the variational symbol. 

Let us now assume that a is function of a real parameter y; -» < y < +«, 
in addition to x^, x^, and x^; a s a(x 1 ,x 2 ,x 3 ; y). Moreover, we consider a 

-* 

proportional loading PPj(x), where p is a finite, dimensionless, real number; 

0 £ p < ®. In this way, the plane of p - y is divided into regions of stabi- 
lity and instability by equation (5.13). The effect of the linear viscosity 
(Eq. 5.3), in this case, is to make the stability regions a closed set (ex- 
cept, possibly, for a set with measure zero; a finite number of isolated 
points in this plane). 

The limiting condition for the flutter of system (5.1), by contrast, is 
obtained when 

2tt 

(i) 

S * J o J v tc ijkje“i,A,i - dv dt = 0 


where m is the frequency of steady state oscillation of the solid about its 
unperturbed state. The motion of the solid decays if H. > 0 and amplifies 
if H 3 < 0. J 

For continuous systems with slight damping, Nemat-Nasser £53] proved 
that the flutter load parameter of the undamped system is an upper bound for 
that of the system with slight damping. He also established the necessary 
and sufficient condition for stability with respect to flutter and suffici- 
ency conditions for stability with respect to divergence and flutter. Based 
on energy considerations he further suggested in [53] an approximate method 
of stability analysis which reduces to the usual energy criterion for the 
case of conservative loading. A complex treatment of a class of one-dimen- 
sional continuous systems was suggested in [54]. 

For a further discussion of stability of continuous (not necessarily 
nonconservative) systems, reference should be made to the recent work by 
Hsu [55] and by Knops and Wilkes [56], 
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CHAPTER VI 


METHODS OF ANALYSIS 


6. 1 Discrete Systems 

6.1.1 Introduction 


The mathematical analysis of stability of discrete systems is most readily 
carried out using Routh-Hurwitz criteria to determine the nature of eigenvalues, 
a subject amply explored in the literature, see, e.g., [18]. As an alternative, 
one could think of applying Liapunov's direct method, [14], with suitable modi- 
fications. In conservative systems this method is tantamount to the well-known 
energy method, but in systems with follower forces suitable generalizations are 
required. Such a generalization has been presented by Walker [57]. The ap- 
proach has the advantage, over an eigenvalue analysis, that the manner, in which 
various parameter changes influence stability, becomes much more visible. Below 
the idea and examples presented in Ref. [57] are reproduced. 

Let us examine first a conservative discrete dynamic system in vector form 

tfq + Cq + Kq = 0 (6.1) 

l~t/ —t* 

where q is an n-vector of displacement, M, C, K are n x n matrices, and M is 
positive definite and symmetric. Most early work is based on energy conside- 
rations, and considers the so-called "conservative" problem (K symmetric and 
positive definite) with "dissipation" forces (C symmetric and positive semi- 
definite) and/or "gyroscopic" forces (C skew symmetric) [58]. For this problem, 
the total energy 

E = ~ q T Mq + ^ q T Kq (6.2) 

is a positive-definite function, having the time derivative 

E - - q T Cq (6.3) 

Depending on whether C is definite, semidef inite , or zero, various exact 
statements can be made concerning stability or asymptotic stability of the 
equilibrium q » q * 0 [58]. 

6.1.2 A "Generalized Energy" Function 

The generalization in Ref. [57] starts by defining n-vectors q^ 

^2 * q and placing the system in the form 

ii -& 

i* ' - x' 1 *3.1 - M ' 1 C S2 


■ q> 

(6.4) 
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while assuming M =» M^, det |m| ^ 0, det |k| ^ 0. These assumptions are main- 
tained throughout. 

Consider now the general quadratic form 


v - Si* F 3i + it G iU + & "Si < 6 - 5 > 

T T 

where F, G, H, are n x n matrices, F * F , G * G . Taking the time derivative 
according to the equations of state, we have 

V = - qj kV 1 Hq x - q^ [2GM _1 C - lfj ^ (6.6) 

+ q 2 X [2F - 2GM -1 K - C T H _1 hJ q^ 


The above function and its derivative are sufficiently general for a full 
attack on the stability problem by the methods of Liapunov, but are also too 
complicated to be of much practical value. It seems desirable to sacrifice 
some generality if a significant reduction in complexity would result. Keeping 
in mind that when the energy method works, it works very easily, suppose we now 
restrict the function V to have two of the properties which the energy function, 
when it is meaningful, normally has; namely: (1) If C 3 0, then V = 0, (2) If 

C 3* 0, then V depends only upon q^. These conditions imply 


F = GM _1 K = [CM" 1 K~J T 
H = 0 


(6.7) 


Thus, the restricted function is 

V = q/ GM _1 K q L + q/ Gq^ 

where G and GM are symmetric, and its derivative is 


( 6 . 8 ) 


V =• - 2 ^ GM _1 C ^ (6.9) 

He note that GH ^ C need not be symmetric, and that ^ s 0 if GH is skew- 

symmetric. We also note that for symmetric K we may choose G “ M and produce 

the energy function, although we need not do so. Since the symmetry require - 

-12 2 
ments on G and GM K result in (n -n) linear equations in the n elements of 

G, there are normally n independent matrices G (producing n independent func- 
tions V) which satisfy these requirements. 

We are now in a position to draw some conclusions: 
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Theorem 1 . If there exists a matrix G such that GM -1 C is positive de- 
finite, while G and GM are symmetric and positive definite, the equilibrium 

is asymptotically stable* 

Theorem 2 . If there exists a matrix G such that GM^C is positive semi- 

definite, while G and GM are symmetric and positive definite, the equili- 
brium Is stable. 

Theorem 3 . If there cueists a matrix G such that GM -1 C is positive de- 
finite, while G and GM -1 K are symmetric but not both positive semidef inite ,* 
the equilibrium is unstable. 

Theorem 4 . If there exists a matrix G such that GM *C is skew- symmetric ,** 

while G and GM are symmetric and definite of the same sign, the equilibrium 
is stable but not asymptotically stable. 

Theorem 5 . If there exists a matrix G such that GM~*C is skew-symmetric, 

-1 T -1 T 

while G and GM K are symmetric, the function V * GM K is an 

integral. 

Although the conditions of Theorems 1-4 are only sufficient, not neces- 
sary, it seems that in the great majority of problems one of the first three 
theorems should prove applicable. A definiteness requirement on an n x n 
matrix results in n inequalities and, as previously noted, symmetry require- 

-12 2 
ments on G and GM K result in only (n -n) equations in the n elements of G. 

Theorems 4 and 5 are les6 generally applicable, particularly when C ^ 0, since 

1 2 

their satisfaction may involve up to (3n -n) equations and, for Theorem 4, 

up to 2n inequalities. However, it may be noted that every result previously 
obtained by the energy method [53](K is then necessarily symmetric) is included 
here by letting G M. 

There are normally n independent matrices G which satisfy the symmetry re- 
quirements on G and GM *K. If the problem is specified in terms of parameters 
and is being attacked analytically, there is often an apparent choice for the 

n arbitrary elements of G which simplify the definiteness conditions on GM 

or, if GM is skew-symmetric, the definiteness conditions on G and GM 
When M is diagonal, a usable result is sometimes found by setting diag [G] * 
diag [M], thus specifying n elements of G a priori. This choice is one which 
produces G C H in the absence of follower forces, and this may be desirable 
since the energy function usually works well in such problems. 


A definite matrix is also semidef inite, as is the zero matrix. 
The zero matrix is also skew-symmetric. 
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Two examples Involving follower forces are presented to illustrate the 
use of these procedures. 


Example 1 , Consider the system 

M i + C 1 + K SL = 0 


where 


M 


0 

0 m 0 


0 0 
0 0 


If we choose diag [G] = diag [M], G « G , 


G = 





( 6 . 10 ) 


( 6 . 11 ) 


( 6 . 12 ) 


and 


GM _1 K = 


kl+ k 3 g /® 2 

k lS /m l +k 3 


-k 3 g /“ 1 +k 2 


For GM K to be symmetric, we 


find 


(6.13) 


" k 3 + k 2 8 ^ m 2 * k l 8 ^ m l + k 3 

(6.14) 

which implies 


2k 3 m l m 2 

8 ‘ k 2 m r k i m 2 

(6.15) 

Since GM ~ 0, we will utilize Theorem 4. Conditions 

definiteness of G and GM are, assuming >0, > 0, 

for positive 

2 

(a) - g >0 

(6.16) 

(b) + m 1 k 2 > 0 

(6.17) 

(c) (k 3 2 + g 2 ) >0 

(6.18) 
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By Theorem 4 these are sufficient conditions for stability, but they also hap- 
pen to be the necessary and sufficient conditions for distinct purely imaginary 
eigenvalues. 

A result which is immediately apparent from this analysis is that propor- 
tional damping always results in asymptotic stability when the preceding con- 
ditions are satisfied. That is, if 

C « + 3K(<*p * 0, a + p > 0) (6.19) 


then 

GM _ 1 c - qG + p GM -1 K (6.20) 

and GM~^C is positive definite since G and GM~*K are positive definite. 
Theorem 1 thus implies asymptotic stability. 

This is not to imply that C need be proportional, or even dissipative, 
for an answer to be obtained. Consider the general matrix 


C ■ 


c i 

C 3 



for which, using the previous G, 


GM -1 C 


’ C l + c 3 S /m 2 

c 3 + Cjg/nij 


we have 
c 4 + c 2 8 /m 2 
c 2 + 


( 6 . 21 ) 


( 6 . 22 ) 


Assuming conditions (a), (b) , (c) , are satisfied, Theorem 1 implies asympto- 
tic stability provided 

(d) c i m 2 + c 3& > 0 (6.23) 

(e) (CjC 2 - c 3 c 4 > * ( m i m 2 - 8 2 ) > 0 (6.24) 

If (d) and (e) are satisfied, but one or more of (a), (b), (c), is 
strictly violated, Theorem 3 implies instability. Stability is assured by 
Theorem 2 when (a), (b), (c) are satisfied while (d) and (e) are weakly sa- 
tisfied. 


For certain values of the parameters, it is clear that a definitive re- 
sult is not given by the preceding analysis. However, all of the preceding 
was based on one specific choice for diag [G]. Other choices can be made and 
results obtained which are applicable under different parameter restrictions. 

It may be noted in general that for problems having C ■= 0, Theorem 4 is 
the only one of the first four theorems having a possibility of success, and 
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it can be satisfied only for systems which would be asymptotically stable with 
the addition of proportional damping. 


Example 2 . Let us consider the problem discussed in Sect. 3,1 which should 
serve as a fair demonstration of the operation of the proposed method. The 
equations of motion are 


l4q + Cq + Kq - 0 

ro 


(6.25) 


where 


“3 1~ 


"0 0 ” 


2-y ory-l 


, C = 


, K = k 


_1 1_ 


_0 °_ 


-1 1- (l-ar)Y 


(6.26) 


2 

where m£ > 0, k > 0, and y =■ Pi/k. The parameters of interest are y and as 
which relate to the magnitude and direction of the nonconservative load. 


Choosing a general symmetric G, 


cf l g 


g 


we find that the symmetry of GM K implies 

g [1 - 2 (1 -cOy] - <* 2 (Y-5) - flf^Y-2) 


Since GM C 3 0, we hope to apply Theorem 4. 


(6.27) 


(6.28) 


Looking ahead to possible algebraic difficulties in the definiteness con- 
ditions on G and GM *K, let us try the simplest choice for (c^, , i.e., one 
which produces g - 0: 


Thus, 


and 


* 5 - y 
= 2 - y 


G 


5-y 0 

0 2-y 


GM _1 K 



(5-Y) (3-Y) 

- (5-Y) (2-Y) 


-(5-y)(2-Y) 

(2- Y )[4-(3-2cr)Y] 


(6.29) 


(6.30) 


(6.31) 
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We see that G is positive definite for Y < 2. Assuming this restriction, we 
find that GM is positive definite provided 


1 + (l-or) (y 2 - 3y) > 0 


(6.32) 


The bifurcation plot in the (Y,ar) plane is shown in Fig. 6.1. We have just 
obtained the shaded region below y * 2 as a region of stability by Theorem 4. 
However, eigenvalue analysis [22] shows the entire shaded region to be the 
exact region of stability. Since our theorems are only sufficient, not neces- 
sary, it may be worth-while to try another choice for G. We see that the pre- 
vious choice for G fails at y * 2 because m 2 “ V is required to be posi- 
tive, a necessary condition for positive definiteness of our previous G. Let 
us simply reverse our definition of end choose 


flf 1 ■ 5 - Y 
* 2 - Y - 2 


(6.33) 


for which 


g - 2(Y-2)(Y-5)/Ll - 2<1 -c0y] (6.34) 

Some tedious calculations show that the conditions for positive definiteness 
of G and GM "St are now 


2 < Y < 5 

41 - 32y + 4ory (1+ory) + 8 y 2 (1-<*) > 0 
7 - 4y + 2 ory > 0 


(6.35) 


1 + (l-or) (Y -3 y) > 0 

These conditions define the shaded area above y = 2 in Fig. 6.1. Thus, with 
the exception of the line Y ® 2 , we have determined the entire region of the 
parameter plane which produces a stable equilibrium. 

We may perform a special investigation for Y * 2. We have then 


G - 


— 1 

Q 

OQ 

1 

9 

gm-'k - “ 2 

Ofj- 3g g(4or-2) 

L 8 


4 

g-3a 2 a 2 (4or-2) 


(6.36) 


r 1- 


where the symmetry of GM K implies 


3q? 2 

3-4or 


(6.37) 


95 



Choosing of- =* 1, we find that the conditions for positive definiteness of G 
-1 1 

and GM K are now that there exist an > 0 such that 
(3-4or) 2 - 9a 2 > 0 

(4o'-2)[(3-4of) 2 - 9 c* 2 ] > 0 (6.38) 

(3-4a) 2 - 90*2 (3 -4a) > 0 

2 13 

These can be satisfied by choosing 9 ot^ < (3-4o0 , provided a > j, t This 

completes the determination of the entire region of stability in the parameter 
plane. Again we note that the addition of any form of proportional damping 

C * TIM + pKOlg £ 0, Tj + p > 0) (6.39) 

leads to asymptotic stability of the equilibrium by Theorem 1, provided (y,a) 
is in the shaded region of Fig. 6.1. The difficulty of reaching this general 
conclusion by eigenvalue analysis need not be dwelt upon. 

6.1.3 A General Approach 

When problems are encountered for which there seems no apparent choice 
for the n arbitrary elements of G which produces a usable result, or in which 
the order is so great as to require the use of a computer, it may be desir- 
able to apply the method in its fullest generality. A systematic approach 
might be as follows: 

1) Define diag [G^ 3 4 5 , where e^ = (1,0,...), “ (0,1,0,...) etc. 

2) Determine G^ by the symmetry of G^ and G^M 

n 

3) Set G = ^ a iPi 

1-1 n 

4) Note that GM _1 K o^GjM^K 

i=l 

n 

GM _1 C = ^ cr i G i M -1 C 
i=l 


5) Determine whether a vector £ exists such that GM is positive defi- 
nite or semidefinite , if so, use Theorems 1, 2, or 3. 
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6) Determine whether a vector cy exists such that GM is skew- symmetric ; 
if so, use Theorem 4 and/or~5. 

6.1.4 Exceptional Cases 

Although it is true that there are problems for which there exists no 
matrix G satisfying the conditions of any of the first four theorems, implying 
stability cannot be determined in this manner, it seems that there is often 
something quite unusual about such problems. One such is the following: 

Hq + Cq + Kq ® 0 

^ ^ rw 


M - 


where > 0, > 0, and f k 2 m l* 

Using the most general form for symmetric G, we have 


m^ 0 


m 0 


0 c 
-c 0 


0 

k„ 


(6.40) 


and 


“i 8 

g O'- 


GM _1 K 


m l m 2 


Symmetry implies 
g = 0 

and therefore 


a l k l m 2 


gk 2 m 1 

a 2 k 2 m l 


GM _1 C = c 


m l m 2 


0 

-a„m 


2 1 


a l m 2 


(6.41) 


(6.42) 


(6.43) 


(6.44) 


and it is impossible to satisfy Theorems 1 or 3 for any choice of (a^, a^). 

Theorem 4 can be satisfied by choosing or^ ■ m^, a ^ provided GM 

is then positive definite. Since 



fm 0 “I 


, 

k i 

0 “ 


1 

o 

1 

& 

» 

GM K = 

0 

k 2 


(6.45) 
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we see that Theorem 4 implies stability provided >0, > 0. 

eigenvalue analysis assures stability provided only that 

+ k 2 m 2 > 0 
k l k 2 > " ° 2 


However, 


(6.46) 


which allow the possibility of one of the k. actually being negative. However, 
-1 1 

noting that our GM * I, we see that Theorem 3 assures that instability will 
occur if C is increased by any positive definite matrix and either of the k^ 

is negative. That is, when the k^ are in the range for stability predicted 

by eigenvalue analysis, but not in the range for which any of our stability 
theorems apply, the addition of any complete dissipation leads to instability. 
On the other hand, when the k^ are in our allowable range, k^ > 0, k 2 > 0, the 

addition of any complete dissipation leads to asymptotic stability by Theorem 
1, a much more natural result. 

6.1.5 Remarks 


Although each of the examples was begun under the assumption of no dis- 
sipation and we initially applied Theorem 4, this was only for the purpose of 
illustrating that we can make qualitative statements about the effects of 
various types of damping via Theorems 1 and 3. In practice. Theorem 4 is the 
most difficult to use and, in the real world, it is the least likely to apply. 
When the linear approximation indicates stability which is not asymptotic, 
the stability or instability of the physical system is determined by parameter 
errors and/or any slight nonlinearities [58,15]. Therefore, Theorem 4 is use- 
ful in problems which are primarily of academic interest. 

These remarks do not generally apply to Theorem 2, although it too con- 
cludes only stability. If the conditions of Theorem 2 are satisfied and those 
of Theorem 1 are not, one can often still conclude asymptotic stability by use 
of an invariance principle due to LaSalle [59]. This usually requires a de- 
tailed analysis of the differential equations, however.** 

Theorem 5 allows the rapid generation of up to n integrals, particularly 
when C = 0, and so permits a possible reduction in the order of the system 


This is a well known result concerning the phenomenon of gyroscopic 
stabilization [58]. 

In particular, one must determine whether any invariant set other than 

T -1 

( q l, q^ = (0, 0), is contained in the set defined by GM C %£ * If 

not, the conditions of Theorem 2 imply asymptotic stability of the equilibrium 
[59]. 
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from 2n to as low as n. However, the algebraic difficulties of such reduction 
may be considerable* The primary use of Theorem 5 is to determine "constants 
of the motion," in problems in which these are of some interest. 


6.2 Continuous Systems 
6.2.1 Introduction 


The class of nonself -ad joint eigenvalue problems is not as extensively 
investigated as that of self-adjoint ones and special attention must be paid 
to the mathematical methods used in determining the eigenvalues. A justi- 
fication of applying the Galerkin method to bars subjected to nonconservative 
loads was offered by Leipholz [60] and the convergence of this method, as ap- 
plied to the same problems, was studied by him in Ref. [61]. Leipholz also 
used the method of finite differences in evaluating the eigenvalue of an 
elastic bar subjected to a uniformly distributed tangential load [62]. He 
extended the applicability of the Galerkin method to a broader class of li- 
near nonself-adjoint eigenvalue problems [63] than those studied in Ref. [62], 
This class contains all eigenvalue problems which arise from originally self- 
adjoint problems by addition of a linear differential expression which de- 
stroys the former self-adjointness. The convergence of the Galerkin method 
for nonconservative stability problems of plates and shells was studied by 
Leipholz in Ref. [64]. A discussion of the Galerkin method as applied to 
systems with damping is discussed by Leipholz in Ref. [65]. 

Levinson [66] has shown that for certain problems the Galerkin method 
converges for a broader class of trial functions than assumed by Leipholz. 
Further, Levinson extended Hamilton's principle and the Ritz method such as 
to make them applicable to nonconservative problems. 

The Ritz method in nonvariational formulation was applied to nonconser- 
vative problems by Marchenko [67]. Both the Ritz and the Galerkin methods 
have been extended further by Leipholz [68]. In particular, the conditions 
to be fulfilled by the coordinate functions are weakened; these functions 
need not satisfy the dynamical boundary conditions and under certain circum- 
stances not even the geometrical ones. This study includes also some consi- 
derations of convergence. 

In treating dissipative dynamic systems of mathematical physics, which 
are governed by nonself-ad joint linear operators, it is often found conveni- 
ent to introduce the adjoint system (or field) and to consider formally a con- 
servative process [69], [70]. The original field contains an energy sink, and 
in the adjoint field an energy source of the same strength is incorporated in 
order to make the combined field conservative. 

It is of interest to note that the notion of the adjoint field can be in- 
troduced also in treating nondissipative, nonconservative systems, i.e., dy- 
namic systems subjected to circulatory forces. In particular, in structural 
systems subjected to follower forces, the consideration of adjoint force fields 
leads to interesting consequences. Indeed, for this class of nonconservative 
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systems, both the original field and its adjoint force field are associated 
with energy sources [71], and yet the combination of these two fields results 
in a conservative one. 


As an example, consider the Beck problem [33], i.e., a cantilevered elas- 
tic bar subjected at its free end to a compressive follower force (see Fig. 
6.2). The equation of motion and the boundary conditions are 


.4 .2 ^2 

+ F ZJL + *JL 

4 2 2 

Bx 4 dx Z dt^ 


0 


7 .g.0 


at x = 0 


(6.47) 


&-&-0 at x = 1 
dx dx 


where dimensionless quantities are employed. We now construct the adjoint 
boundary-value problem by considering a function z = z(x,t), defined for 
0 £ x £ 1 and t ^ 0, such that the following equation of motion and boundary 
conditions at x = 0 are satisfied identically: 


^4 ^2 .2 

a z . v a_z a z 

ax 4 ax z st* 


o 


z 



at x = 0 


We then seek boundary conditions for z, at x B 1, such that 



(6.48) 


(6.49) 


If we now integrate the right side of the preceding equation by parts and use 
boundary conditions in Eq. (6.47), we immediately obtain the following boun- 
dary conditions for z at x = 1: 


2 3 

+ Fz = 0 ^-f + F |^ = 0 at x = 1 (6.50) 

dx 2 fcc 9X 


Equations (6.48) and (6.50) now define the system adjoint to the Beck problem. 
The inspection of Eqs. (6.48) and (6.50) readily reveals that they describe 
the Reut problem [8] sketched in Fig. 6.3. 
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It was shown in Refs. [8] and [33] that both systems depicted in Figs. 

6.2 and 6.3 lose stability for the same value of the load, i.e., F fir * 20.03. 

For F > F cr , energy is transferred to the bar by the work of the applied 

force, which in turn increases the conservative energy of the system, making 
the response unbounded (flutter). Thus, both force fields are associated with 
source of energy. However, it can be readily seen, both on mathematical and 
physical grounds (see Fig. 6.4), that the combined system is conservative, 
and the bar shown in Fig. 6.4 is incapable of losing stability by flutter. 

Indeed, the loss of stability in the combined case will occur by divergence 

(buckling, attainment of another equilibrium configuration). In conclusion, 
it should be mentioned that adjoint systems can be constructed also for the 
nonconservative problems discussed in Refs. [8], [72] and [73]. 

Adjoint systems have been also examined in Ref. [74]. The possibility 
of constructing adjoint equations for the purpose of developing approximate 
methods in aeroelasticity similar to energy methpds was indicated already in 
Ref. [75]. The usefulness of adjoint systems in solving stability problems 
of elastic continua with follower forces was exhibited in Ref. [76], as de- 
scribed below. 


6.2.2 Stability of an Elastic Continuum 


Let us consider an isotropic, homogeneous, elastic solid occupying a vo- 
lume V bounded by a finite surface S. It will be assumed that on one part of 
the boundary of the solid the displacements are prescribed so as to pre- 
clude a rigid body motion. The body is at rest and in a state of initial 
stress cr^ j , i,j * 1,2,3, due to the applied nonconservative (follower) forces 


on the surface S - S q of the solid. To study the stability of this rest po- 
sition the system is slightly perturbed and the type of ensuing motion is 


Referred to an orthogonal cartesian coordinate system x^, Bolotin 


studied. 

[8] has obtained the following equations for the ensuing motion: 

A, 

0 in V 


_JL f|k \ + R _B_ / *l \ 

dXj V^ijia Bx^ / P dXj V jk Bx fc / p ^ t 2 


(6.51) 


K l 


50 , 


X ijk/ Bx, n j + Pa jk Bx k n j 


Pp. on S - S 
i o 


(6.52) 


u. B 0 on S 
1 o 

(6.53) 

X ijkJ5 “ X6 ij 6 kX + 2M,6 ik 6 j/ 

(6.54) 

fo. It 3 

‘ k i., 
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In Eqs. (6.51)-(6.54) , p is the mass density, is the displacement 

vector measured from the undisturbed state and n^ is the outward positive 

unit normal vector to S. No body forces are assumed to be present and 3 is 
a parameter associated with the magnitude of externally applied surface trac- 
tions* In Eq. (6.54), X and ^ are Lamp's constants of elasticity. The re- 
peated indices are summed over the range of their definitions and p^ are the 

components of perturbations of the applied surface tractions and their forms 
will depend on the behavior of the nonconservative forces. They will gene- 
rally be homogeneous functions of displacements and their derivatives with 
respect to both space and time. In the present study, however, it suffices 
to restrict to the following expression: 


3u i 

a, .u. + b. -r — 
ij J J Sx 


on S - S 


(6.55) 


where a. . and b, are coefficients which are independent of the vector u. and 
iJ j J 

its derivatives but in general are functions of spatial coordinates x^ . 

We may assume a solution of the above boundary value problem in the form 


icut 


u j (Xi,x 2> x 3 ,t) * u j (x 1 ,x 2 ,x 3 )e 
which results in the following eigenvalue problem: 


(-D 


1/2 


( x ukx ^ ) + s ^ (“jk ^ ) - Au i ■ 0 i * T 


(6.56) 


ijkX bx. 


sT “j + ^jk aT ”j ' e ( a ij u j + ‘it) °" s - s. «•”> 


u^ = 0 on S q (6.58) 

A = - CD 2 (6.59) 

Equations (6.56) -(6.58) constitute a nonself -ad joint homogeneous system 
and stability of the solid will be governed by the character of the eigen- 
values A m , m * 1,2,...», for nontrivial solutions. In view of the fact that 
the applied surface tractions are not derivable from a potential, it is not 

possible to express the eigenvalues A™ in the form of a ratio of two positive- 
definite integrals, and thus the usefulness of variational principles seems 
dubious in this case. 
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6*2.3 The Adjoint System 


By constructing an adjoint system by means of certain mathematical rela- 
tions analogous to the definitions in the theory of ordinary differential 
equations, A may be expressed in terms of the original and the adjoint vari- 
ables, and as a consequence A will assume a stationary value. In the theory 
of ordinary differential equations, a system adjoint to one governed by a dif- 
ferential equation and boundary conditions may be constructed formally by re- 
peated integration by parts [77]. Being guided by this observation we examine 
the problem 




^ * 

\ 

ijkX 3x, ) 


p ( a jk ^ 


0 in V 


du* 

X ijk£ 3^ n j + 3^ n j “ P a ij u 5 + P c ij u 5 on S - S , 

u* • 0 on S 
i o 


(6.60) 

(6.61) 

(6.62) 


as being possibly adjoint to that given by Eqs. (6 .56) -(6.58) . Here, c^ is 
a function of bj , u^ and its derivatives. If an adjoint system is to be de- 
fined through equations (6.60) -(6.62) , one must obtain c ^ by solving a cer- 
tain homogeneous integral equation on the surface S - S q . The above-mentioned 
integral equation reduces to satisfying the following: 


j dx. 


U j C ji 


(6.63) 


Expression (6.63) involves three independent equations in nine unknown 
quantities c^ and thus an adjoint system is not uniquely defined [77]. As 

a consequence of Eq. (6.63) the following holds: 



(6.64) 


This expression appears to be similar to Maxwell's reciprocity relations 
in conservative systems, in which case u^ = u*. The bracketed terms are re- 
cognized to be resultant forces associated with the original and the adjoint 
systems, respectively. 
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be the eigenvalues of equations (6.56) -(6,58) , 


Now let A , m * 1,2,...®, 
and A*®, m *= 1,2,...®, those of equations (6.60) -(6.62) , while the correspond- 
ing eigenfunctions are u™ and u* m , respectively. From Eqs. (6.56), (6.60) 
and (6.64), we have 


■ 4 

[^< 

^ m 

' Su k ' 

, X ijkje bx £ j 

) + B 

= r u m 

PL ( 

' K n > 

) + B ^I 

j, u ‘ 


w ijk£ > 

-A* n r 

Jtt 

m *n 

u u i dv 



m 


jk 


du 5 




dV 


(6.65) 


Therefore, 


, A m A *n x P m *n . 

(A - A ) J u ± u £ dV = 0 


( 6 . 66 ) 


At this point we wish to apply the argument of Roberts [77] to prove 
that the sets of eigenvalues {A m } and {A* m } are identical. Let us suppose 
that {A m } and [A* n ] are not identical sets, then 


J v <$v» - •» J v »r (* 


ijk/ dx 


du 


m 




dV = 0 

(6.67a) 


and for the special case when m = n. 


f u*V m dV = 0 (6.67b) 

If the set of eigenvectors [u^ m ] is complete, Eq. (6.67b), together with 

jLm 

Eq. (6.67a), would imply that u^ is identically zero, which is not nontrivial. 

Hence the two sets of eigenvalues are identical. Also, similarly to the pro- 
perty of orthogonality of principal modes in the theory of small vibrations, 

Eq. (6.67a) reveals that the two sets of eigenfunctions [u^ m ] and {u* m } are 

bi-orthonormal, i.e., each function of either set is orthogonal to every mem- 
ber of the other set except those which belong to the same eigenvalue. 
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From (6.65) it also follows that 


m I v “l [ok. ( X ijki bx t ) + e Ok. ( 5 Jk 0x fc )] dV I ; 

4 K m “*r « 


(say) 


( 6 . 68 ) 


.HI 


Let us consider now the effect on A due to infinitesimal variations 
and 6u* m which are arbitrary except that they satisfy the boundary con- 
ditions (6.57), (6.58) and (6.61), (6.62). Therefore, 


6A m - — (61 A m 6I 2 ) 


dV “V 


L { ( x m 


V 

. ijki dx. 

J ** 


du. 




Am 


r m 1 

u. u 

♦ «r 

A m ( m c Am Am- m\\ 

- A l u i 6u^ + u ± JY dV 


(6.69) 


Equation (6.69) reduces, after application of the divergence theorem 
and satisfaction of boundary conditions, to 

J v { -r (‘ 1)M ) + ' ^ ^ ¥ )• >V] 


LA- 

*m 

(. - 

\ 4- Q & 

/ i N . .« *mTl 

Lax. 

\ ijia sx. , 

) + P dx. 

V jk dx. ) U i Jj 


dV 


(6.70) 


Equation (6.70) is clearly a useful version of a variational principle 

m 


and implies that if Eqs. (6.56) and (6.60) are obeyed, 6A is zero with an 

m Am 

accuracy of first order for all small arbitrary variations 6^ and 6u^ 
that satisfy the boundary conditions (6.57), (6.58) and (6.61), (6.62), 
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respectively. Thus a definite statement can be made regarding the error in- 
volved in stipulating that the eigenvalues are stationary values. 

6.2.4 An Approximate Method of Stability Analysis 

The extremum property of the eigenvalues A m , as expressed by Eq. (6.70), 
suggests an approximate procedure for their determination, in the spirit of 
approximate methods for self-adjoint systems based on variational principles. 

We may select two sets of trial functions ^”(0^, a^, • • •) and 

which satisfy the appropriate boundary conditions and contain undetermined 

TO 

parameters and oty An approximate expression of the eigenvalues A is 
obtained, by using Eq. (6.68), as a function of these parameters. A station- 
ary value of A™ is then obtained by determining the parameters from equations 
of the type 


3 “j 


0 ; 


K 


which is reminiscent of the Rayleigh-Ritz procedure for conservative systems. 

6.2.5 Illustrative Example 

In this section we wish to apply the approximate method discussed above 
to investigate the stability of equilibrium of a cantilevered bar subjected 
to a follower load. The governing equations of motion may be expressed as 
[ 8 ] 


.4 .2 

d u , _ d u 

4 + r 2 

dx dx 


2 

u> u 


0 ; 


0 £ X S 1 


(6.71) 


T* - 0 

dx 


d 2 u 

dx 2 


A 

dx 3 


at x = 0 


at x 


(6.72) 


In Eqs. (6.71) and (6.72), dimensionless quantities are employed and cu 
denotes the frequency of oscillation. The equations of an adjoint system of 
this problem, which was first discussed in Ref. [78], are as follows: 


A * a 2 * 0 

d u . „ d u 2 * 
7— + F U) U 


dx 


dx 


(6.73) 
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* 

u 



at x » 0 


d it* 
dx 2 


+ Fu* 


du 

dx 3 


+ F 


dx 


at x 


(6.74) 


2 

The eigenvalue oo in the two problems will be the same as established in 
general in the previous section, and we wish to determine it approximately. 

Ve assume, then, that u and u* may be written in the form: 


N 


y or u 

L n n 


n*l 

N 


* V * * 

tf = ) O' U 

/ i n n 


n-1 


(6.75) 


(6.76) 


where u^, u* are certain assumed functions of x which satisfy the boundary 

conditions (6.72) and (6.74), respectively, and ar n , o'* are constants to be 

determined as discussed. We multiply (6.71) by u* and integrate over the 
length. If we substitute the expansions (6.75) and (6.76), the following 
relation is obtained: 


(A) 


1 N 4 2 N * 

I y o*u* (L + F £--) faudx y *** A 

J L m m \, 4 , 2 / L n n L m n mn 

° m=l dx n«=l m.n°l 


where 


(6.77) 


r y c^* u* y » u dx 

J L m m L n n 


m=l n=l 


1 

m,n=l 


B 

n mn 


A A 2 

A u d u 


,u u u u v 

A ■= f u* ( T” + F — 2 s ) dx 
mn o ” v dx 4 dx 2 ' 


mn 


f u*u dx 
J m n 


To obtain the best possible result, we must now seek an extremum of the 
2 

expression for ou considered as a function of the parameters and c? n . A 
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2 

simple and familiar way would be to treat u) as a Lag rang i an undetermined 
multiplier and seek directly the stationary value of the following: 


I 


N N 

"r* V A - y cPa B 

2 Zj m n mn L> mn mn 

^ m, n=l m,n“l 


(6.78) 


by requiring that 


31 s ai 

Bo* dc* 
m n 


0 


Since u and u* are functions that satisfy the adjoint relations in the 

sense discussed before, it is a simple matter to show that dl/do?* and dl/do' 

mi 

would result in two matrix relations which are adjoint to each other and 
thus they would yield identical eigenvalues. Therefore, in the sequel only 
the following relation will be considered: 


si 

So* 

m 


= 0 


(6.79) 


Equation (6.79) is a homogeneous, linear, algebraic equation in and, 

therefore, a nontrivial solution exists only if the determinant formed by 

the coefficients of a vanishes. This results in a polynomial equation for 
2 n 

m which represents approximately the frequency equation of the system. 

Let us consider the following specific trial functions with N = 2: 

/ 2 2 3 . x 4 \ , / 3 4 3 5\ Qns 

u = Qfj ( x - 2 X + ■g- 1 + 0f 2 ( x - X + jp x ) (6.80) 


»* - rf { x 2 - + -riu2_ } 

L F 2 + 6F + 72 F Z + 6F + 72 J 


+ * f 3 . 2(F 2 +12F + 120> x 4 + 

^ -w ^ 1 _I_ "> /. A 


F + 16F + 240 


F + 6F + 72 
F 2 + 16F + 240 


} 


(6.81) 


Functions (6.80) and (6.81) satisfy the boundary conditions (6.72) and (6.74) 
respectively. Following the procedure as discussed before, we obtain the 
frequency equation: 
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2(F + 4F + 24) 
F 2 + 6F + 72 

F 2 + 12 
F 2 + 6F + 72 


A z c 2££_± 12F + 120) 
F 2 + 16F + 240 


B' 


F + 6F + 72 
F 2 + 16F + 240 


Equation (6.82) will yield distinct real root s s for vanishing F, and when 
F is increased the two roots will coalesce at the critical value F * Y^ beyond 

which (6.82) will yield complex roots. By trial and error Y^ is computed to 

be 19.45, whereas a more precise calculation by Beck [33] yields Y^ * 20.05. 

Incidentally, if one uses only the trial function (6.80), the method of 
Galerkin yields Y^ *= 20.6. This result was first computed by Levinson [66]. 

A similar approximate method of stability analysis was worked out inde- 
pendently of the above by Ballio [79]. 
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6*3 Energy Considerations 


It appears appropriate to discuss energy considerations in the context 
of methods of analysis, because such considerations, as in the class of con- 
servative systems, may lead to the establishment of approximate methods of 
analysis. In addition, energy considerations may be useful in deriving the 
differential equations of motion (as well as the boundary conditions) and to 
provide additional insight into certain aspects of instability phenomena. 

As compared to the vast amount of literature concerned with stability 
of mechanical systems with follower forces, it is somewhat surprising to 
observe that only few studies contain energy considerations. While investi- 
gating the dynamics of articulated pipes, Benjamin [80] invoked Hamilton's 
principle and discussed the energy transfer to the system. In Ref. [71] an 
extension of the usual energy method was proposed, such as to make it appli- 
cable for the stability analysis of circulatory systems with and without 
velocity-dependent forces. Energy considerations formed the basis of deriving 
equations of motion in systems with follower forces in Ref. [81]. 

Energetic and thermodynamic considerations in stability of conservative 
and nonconservative systems were discussed in Ref. [82]. In Ref. [83] an ap- 
proximate energy method for finding the relationship between the force para- 
meter and the amplitude of steady-state oscillations of nonlinear, nonconser- 
vative, autonomous systems was suggested. Stability criteria on the basis of 
"equivalent energy" conditions were established in Ref. [21]. 
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CHAPTER VII 


POSSIBILITIES OF PHYSICAL REALIZATION 


7.1 Introduction 

It is a peculiar common feature of much published analytical work on the 
dynamics and stability of mechanical systems with follower forces, that the 
possible physical origin of such forces is not mentioned. The follower forces 
are Introduced into the analysis either through a sketch, with forces being 
merely indicated by arrows, or through a specified functional dependence of 
the forces on generalized coordinates. Thus the problem is reduced immediately 
to a mathematical analysis and the relationship to mechanics (as a branch of 
physics or engineering) becomes most tenuous. The motivation for much of this 
type of work appears to have been sheer curiosity in determining the sometimes 
unexpected behavior of an imagined system, rather than an explanation of ob- 
served phenomena. 

This clearly unsatisfactory state of imbalance in the development of an 
area of applied mechanics can be rectified by paying, as a first step, atten- 
tion to the possible physical origin of the follower forces which are intro- 
duced into the analysis and building, as a second step, actual demonstration 
models, to be followed by a quantitative experimental program. 

Let us discuss in this Chapter some possible origins of follower forces. 

If the mechanical system should be able to lose a position of equilibrium 
through oscillations with increasing amplitudes, a source of energy should 
be coupled, through the follower forces, to the system. In one category of 
problems involving rotating shafts this energy is supplied by the driving 
motor and stability is lost by lateral oscillations. This category of prob- 
lems is deliberately not covered in this report. 

In another category, the energy is supplied through a moving fluid to 
the mechanical system. If the fluid surrounds the mechanical system whose 
stability is being studied, the problem belongs to the broad and technically 
most significant area of aeroelasticity. The kinetic energy of a fluid can 
be transferred to the system also through Internal flow in flexible pipes and 
by means of impinging jets. Some of these possibilities will be discussed 
presently. 

It is conceivable that other forms of energy, such as, e.g., chemical 
and electro-magnetic energy, could constitute appropriate sources which, 
under suitable conditions of coupling, could induce flutter- type instabili- 
ties. Among all these possibilities, the author is aware only of some recent 
work on instability (including flutter) of bars induced by radiant heat, as 
mentioned in Sect. 7.4. 
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7.2 Instability Modes of Cantilevered Bars Induced by Fluid Flow Through 
Attached Pipes 


7.2.1 General 

Let us discuss, as an example, the problem stated in the above heading. 
This particular example has been chosen, because various types of instabili- 
ties occur in a richer variety than, e.g., in a single, axi- symmetric flexible 
pipe conducting fluid, a system discussed in Chapter VIII since some actual 
experiments have been reported. It will be shown that a cantilevered bar 
having two axes of symmetry may lose stability by either torsional divergence, 
torsional flutter or transverse flutter, but not transverse divergence. The 
Coriolis forces can have either a stabilizing or a destabilizing effect on 
both the torsional flutter and the transverse flutter, depending upon the 
parameters of the system [84]. Stability of a similar bar subjected to a 
single eccentric follower force was discussed in Ref. [85]. The treatment 
can be considered a special case of no Coriolis forces. 

7.2.2 Derivation of Equation of Motion and Boundary Conditions 

We consider a thin-walled, cantilevered, elastic beam with two pairs of 
flexible pipes, which are attached to the bar at a distance h/2 from the 
z-axis (so that the whole system deforms as a unit) and pump fluid at a con- 
stant velocity U through the pipes, as sketched in Fig. 7.1. We designate 
the length of the system by L, the torsional rigidity by C = GJ, and the 
warping rigidity by C^ *= EC^, [86], and similar to the work of Benjamin [80] 

obtain the equation of torsional motion of the system, using Hamilton's 
principle. With cp(z,t) denoting the angle of rotation at section z and at 
time t, the strain energy of the torsional deformation is [87] 


V 


1 



[CjCcp ") 2 + C(cp') 2 ] dz 


(7.1) 


where primes denote differentiation with respect to z. The kinetic energy 
is 


T 


1 



(7.2) 


where a dot denotes differentiation with respect to time, m is the mass of 
the assembly per unit of length (exclusive of the mass of the fluid), and r 
is the polar radius of gyration of the cross-section of the system. 


The total kinetic energy of the fluid may be obtained by adding to the 
kinetic energy of the fluid contained within the pipes, T^, the change in 

the kinetic energy of the fluid entering and leaving the pipes during a very 
small interval of time At: 
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(7.3) 


T' = T 2 + 2M0 (± U o 2 - \ U^) At 

where T* is the total kinetic energy of the fluid, M the mass density of the 
fluid per unit length of each pair of pipes, U q the outlet velocity vector, 

and U ± the inlet velocity vector. But *= Ui, where i is the unit vector 

in the z -direction, and U q « ? + riUT where t is the unit vector tangent to 

the top (bottom) pipe at z = L, r is the position vector of the top (bottom) 
pipe at z = L, and n is the ratio of the area of each pipe to that of the 
attached nozzle at the free end. Hence, 6T* becomes 

6T ' * 6T 2 + 2MU(r + nUr) • 6r (7.4) 


The components of the absolute velocity of the fluid are y + U(dy/dz) 

in the y-direction, and U [1 - — (y ') ] - w in the z-direction, where w(z,t) 

denotes the average displacement at section z and at time t in the z-direction. 
then becomes (within an additive constant) 

T 2 » 2M J 1 0 y 2 + Uyy Uw) dz 
o 

But y *= (h/2)cp, which yields 

T 2 - 2M J 1 jjg- 9 2 + qxp' - UwJ dz (7.5) 

o 

With j being the unit vector along the y-axis, we have (see Fig. 7.1) 


t = j sin 0 + i cos 6 = j(y') z=L + i 
= cp'(D] j + i 

r = j(y) z=L - i(w) z=L = \j <ptt>] j - [ W (L)] i 

Then 

* h 2 

(r + nUr) • 6? - nU 6w(L) + [q>(L) + nU <p'(L)] 6cp(L) (7.6) 

where w(L) 6w(L) is neglected (being a term of higher order) • The Lagrangian 
now becomes 

L ■= T x + T 2 - W 1 + 2MnU 2 w(L) (7.7) 
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and Hamilton's principle takes on the form 


2 ^ 2 2 

6 f Ldt - F MU [cp(L) + nV cp'CL) ] 6cp(L)dt - 0 
t, **1, 


(7.8) 


where 


W(I») “ J J 1 * r 2 (cp # ) 2 dz 


Carrying out the variations and using integration by parts, we obtain 

S + 1 2 "”” 2 ' 2 - c ] rf + x™ 2 & + (“ 2+ “ r) 7 f ■ • 

OZ oz ot 


* 0; 


z = o 


.2 

if -0 

oz 


<= l ^[“»« 2 ( 2 r 2 -r)-«]S-o 

oz 


z ■ L 


(7.9) 


He now introduce the following dimensionless quantities: 

5 -f 


T ■ t 


L [mr 2 +(h 2 /2)M]L 4 


H 


CL 


“ = r > 




o' 


B 

l+(a 2 /2)p 


, and F = 


MnU 2 r 2 L 2 
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Equations (7*9) then become 



[2F-k] 



a 


d 2 cp 



0 


<P 



at 5 - 0 






5 - 1 


(7.10) 


which are analogous to those obtained in Ref. [72] for cantilevered bars 
subjected at the free end to follower forces, except for the third term in 
the first equation, which is due to the Coriolis acceleration. As we shall 
see in the sequel, this term can have either a destabilizing or a stabilizing 
effect. That is, for sufficiently small Coriolis forces (n large and p' small) 
the system loses stability (by torsional or transverse flutter) under smaller 
F than obtained when n = ® (no Coriolis forces). On the other hand, for p / /n 
sufficiently large, the critical value of F can be increased by increasing 
3 7n. 


We note here that, in torsional instability similar to transverse in- 
stability, the Coriolis forces have an effect similar to that of internal 
viscous damping [32]. That is, although damping (and also Coriolis forces) 
is a dissipating agency, when it is sufficiently small, it may act as a chan- 
nel for the transfer of energy to the system from the source, which is always 
associated with the type of nonconservative forces considered here [71]. 


7.2.3 Stability Analysis 

Frequency equation . We take the solution of system (7.10) as cp(?,T) = 
^(5)e‘ i ‘ UJT and obtain the following eigenvalue problem: 




| + [2F-k] + (iu>) J(P~\ J] a 2 - u> 2 * = 0 


♦ - - 0; at ? - 0 

^-0 

dr 

$♦[»(*- #)-■>]!!-«, 


§ - i 


(7.11) 
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We then let f(§) ■ Ae*^ and obtain 

X 4 ~ (2F-h)\ 2 - 0) ' “) Qf 2 X - u) 2 = 0 (7.12) 

Equation (7.12) is a polynomial of degree four in X and therefore has, 
in general, four complex roots. Let these roots be designated by X,; 
j = 1,2,..., 4. Then, ^ 

*<« ■ l 

1-1 


which may now be substituted into the boundary conditions to yield four homo- 
geneous equations for four constants A.« These equations are 


4 



j=l 


i 


j-i 


Vj 


= 0 


I 


j=l 


2 

vv -° 


(7.13) 


X 


j-1 


a s 


iX. 


= 0 


2 

where T| - F(2-0' /2)- h. System (7.13) has nontrivial solutions if and only 
if the determinant of the coefficients of A.; j = 1,2,... ,4 is identically 
zero, i.e., the frequency equation is J 

2 2 

A = e (X L X 2 +^X 1 X 2 )(X 2 -X 1 )(X 4 -X 3 ) 

i(Xi+Xo) 2 2 

- e (X 1 4t X3^+TlX 1 X 3 )(X 3 -X 1 )(X 4 -X 2 ) 

i(X^+X 4 ) 2 2 

+ e (^i X 4 +%X 4 )(X 4 -X 1 )(X 3 -X 2 ) (7.14) 

cont. 
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i(X 2 +X~) 2 2 

+ c J (x 2 x 3 +tix 2 x 3 )(x 3 -x 2 )(x 4 -x 1 ) 


i(X 2 +X 4 ) 2 2 

- « (x 2 x 4 +m 2 V ( V X 2 )(X 3“ X l ) 


1 ( X 3 +X 4 ) 2 2 

+ e J (X 3 X 4 * +m 3 X 4 )(X 4 -X 3 )(X 2 -X 1 ) - 0 (7.14) 

where X^, X 2> X 3> and \ 4 are defined as functions of to through Eq. (7.12). 

Torsional buckling . To obtain the condition for divergent torsional 
motion, we let co ■ 0 in Eq. (7.12) and obtain X^ 2 * 0, and X 3 4 = ± . 

2 — 2 * 

Then, with k 6rr and F ■ 2F - h ■ yn , Eq. (7.14) reduces to 

a 2 4 V c osn / _ y (7.15) 

(y+6)( 1-cos ttVy) 

which is identical to the equation obtained for the torsional buckling of a 
cantilevered beam subjected at the free end to follower forces [72]. The 
first branch of the torsional buckling, corresponding to the first mode of 
instability, is shown by the solid line in Fig. 7.2. 

— 2 

Torsional flutter . For given <*, 3, n, 6 and F = yn , Eqs_. (7.12) and 
(7.14) yield the frequencies of torsional oscillations. When F is small, 
these frequencies are all located on the left hand side of the imaginary axis 
in the complex iiu plane and the system can perform only damped torsional 
oscillations. 

As we increase F, one of these frequencies approaches the imaginary axis, 
and for a certain value of F, say F^^, Eqs. (7.12) and (7.14) yield a real 

value for to. If we now increase F beyond this critical value, one of the 
roots of (7.14) becomes complex with negative imaginary part. The beam will 
oscillate with an exponentially increasing amplitude. Consequently, we shall 
seek, for given a, 3, n, and 6, values of ao (real) and "F which identically 
satisfy (7.12) and (7.14). This can be done directly with the aid of a com- 
puter. The computer can be instructed to find the roots of Eq. (7.12) for 
specified values of a, 3, n, 5, a), and y> and then 'compute the value of A. 

By varying the value of u; and y systematically, the critical co and y may 
easily be selected which make both real and imaginary parts of A identically 
zero. This is illustrated in Fig. 7.3 where for a * 1.50, 6 * 1.0, 3 ■ 1.0, 
and n » 1, both real and imaginary parts of A = A^ + iA 2 are plotted against 

2 2 4 

the values of co . We see that for y ■ 3.40, and co * 1.13tt , A is identically 

zero. Similar results may be obtained for other values of a, 3» and n. In 

this manner torsional flutter curves may be constructed. The first branch 
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(the only practically significant one) of torsional flutter is shown in Fig, 
7.2 by dashed lines, for 6 * 1, n * 1, and indicated values of g. The solid 
curve for torsional flutter in Fig. 7.2 is the limiting case when n = » and 
corresponds to the torsional flutter of a cantilevered bar subjected at the 
free end to compressive follower forces [72]. 

It must be noted that, even for relatively large values of 0(n**l), the 
Coriolis forces may have a destabilizing effect for certain values of or. 

(For example, for P * 0.1 and 1.0 < a < 1.35, as is seen in Fig. 7.2.) 

Transverse flutter . In addition to torsional buckling and torsional 
flutter, the bar may lose stability also by transverse flutter [87]. The 
equation of motion and the boundary conditions for this case have been de- 
rived by employing Hamilton's principle in [80] and D'Alembert's principle 

with El, cp(z,t) with y(z,t) and 


(2M + m) « 0 

bt* 

y * = 0; at z * 0 

= 0; at z = 0 

bz bz 


in [87]* Here, we may simply identify C- 
write 


4 2 2 

El + 2MriU 2 + 4MU + 






dzdt 


which, by introducing the following dimensionless quantities: 


5 = f > T- t^[- 4 

L VL (2M+m)L 4 - 


Mnl^L 2 


El 


□ _ M , = _2£_ 

P m » P 2 P+1 


reduces to 


+ ^_2 + A n‘i) aea T + -,_2 0 


as a?' 
. a* 


St 


y = 


a§ 


o; 


5 = 0 


e-i 

a ?' 2 
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Equation (7*12) now becomes 


X 4 -2 fl X 2 -»^ Fl )x-» 2 -0 
and equation (7*14) takes on the form 
i(X 1 +X 2 ) 2 2 


A “ e 


x x x 2 (x 2 -x 1 )(x 4 -x 3 ) 


(7.12 ') 


- e 


i(Xj+X 3 ) 2 2 


X^ ^3 


i(X^+X^) 2 2 


i(X 2 +Xo) 2 2 

+ e X 2 X 3 (X 3 -X 2 )(X 4 -X 1 ) 


(7. 14') 


- e 


i(X 2 +X 4 ) 2 2 


x 2 x 4 (x 4 -x 2 )(x 3 -x 1 ) 


i(X 3 +X 4 ) 2 2 

+ e J X 3 X 4 (X 4 -X 3 )(X 2 -X 1 ) - 0 


For a given (3 and n, we now seek values of to and F^ which identically 
satisfy (7.12') and (7.14'). In this manner we obtain the limit of trans- 
verse flutter, as shown by horizontal dashed lines in Fig. 7.2 for 
2 

Elr /C^ * 1.5 and g ** 0.1, 0.2. In this figure, the horizontal solid line 

indicates the limit of transverse flutter for n ** ® [72]. We note that for 
g * 0.5, 1.0, the transverse flutter occurs at y * 12.2, and 15.8 respectively. 
These values are not shown in Fig. 7.2. 

7.2.4 Analysis of Flutter by Indirect Method 

The method used in the previous section for the analysis of flutter was 
a direct one. That is, for a given system we directly obtained the critical 
values of y and w. One may solve the same problem by an indirect method which 
was employed in [87]. 

To this end we let Xj» j “ 1,2,... ,4 denote the roots of Eq. (7.12). 

Then we have 
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x x + x 2 + x 3 + x 4 = o 

X 1 X 2 + X X X 3 + XjX 4 + X 2 X 3 + X 2 X 4 + X 3 X 4 - - (2 F-h) 

XiX 2 X 3 + X x X 2 X 4 + X X X 3 X 4 + X 2 X 3 X 4 “ ^ a a) 

2 

X ^ ^ ^ ^ = * u) 

The first equation in (7.16a) is identically satisfied if we let 
X^ = a - b - c 

X 2 --a + b- c 

X 3 =-a-b + c 

X^ “ a + b + c 


and from the remaining equations we obtain 


2 

a 


+ b 


+ c 


2F-k 

2 


abc = | ,/l 

C»'5)A 

a 4 + b 4 - 

2 2 2 2 2 2 
2a b - 2b c - 2c a 


2 

(U 


We now let 

a = J (p + iq) 

b = | (p - iq) 

and from (7.16c) obtain 

p 2 - q 2 + 2c 2 = 2F - K 

(p 2 + q 2 )c J (\ P 7 a 2 u) 

(p 2 - c 2 )(q 2 + c 2 ) = a) 2 


(7.16a) 


(7.16b) 


(7.16c) 


(7.16d) 


(7. 16e) 
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where p, q, and c are all real. Substituting from (7.16d) into (7.16b) and 
then into the frequency equation (7.14)* we finally arrive at, after a series 
of tedious manipulations, 

A 2 + iA 2 - 0 


where: 

* cp {2 [3p 2 q 2 - q 4 - c 2 (p 2 + q 2 ) + 4c 4 ] + (p 2 - 3q 2 - 4c 2 ) T] } cos p sinh q 
+ cq [2 [p 4 - 3p 2 q 2 - c 2 (p 2 + q 2 ) - 4c 4 ] - (3p 2 - q 2 - 4c 2 ) T| } sin p cosh q 
+ pq {[(p 2+ q 2 )(q 2 - P 2 + 2c 2 )] + (p 2 + q 2 )T| } sin 2c 
A 2 88 £[p q (q - P ) + c (p + q - 6p q ) - 3c (p - q ) - 4c ] (7.17) 

+ [2p 2 q 2 - c 2 (p 2 - q 2 ) + 4c 4 ] T] ] sin P sinh q 

+ pq [2 [-p 2 q 2 + 3c 2 (p 2 - q 2 ) - 7c 4 ] - [p 2 - q 2 - 2c 2 ] T| ] cos p cosh q 

- pq [[p 4 + q 4 + 2c 2 (q 2 - P 2 ) + 2c 4 ] - [p 2 - q 2 - 2c 2 ] T] } cos 2c 

2 

For an assumed value of c and given a and h = 6tt , we may now find p and 
q such that * A 2 “ 0* Then, from equations (7.16e) the corresponding 

values of F, (3, and cu, for a given n, may be computed. 

The above method is an indirect one, as we do not know, in advance, which 
particular problem is being investigated. Moreover, if a computer is being 
used to find values of p and q which satisfy A^ ** A 2 = 0, it is then just as 

easy to employ the direct method outlined in the previous section. However, 

for small values of Coriolis forces, that is for sufficiently small ^(g'/n), 
one may reduce Eqs. (7.17) by neglecting the higher order terms in c and study 
the effect of small Coriolis forces directly. This we shall discuss in the 
following section. 

7.2.5 The Effect of Small Coriolis Forces 

2 

We consider equation (7.17) and by neglecting 0(c) and higher order 
terms obtain 
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(7.18) 


A x - P [2 [3p 2 q 2 - q 4 ] + [p 2 - 3q 2 ] T| } cos p sinh q 
+ q {2 [p 4 - 3p 2 q 2 ] + [q 2 - 3p 2 ] T| } sin p cosh q 
+ 2pq {(q 4 - P 4 ) + (p 2 + q 2 ) 1] } ■ 0 

^2 * pq [(q 2 - p 2 ) + 2T|] sin p sinh q 
2 2 2 2 

+ [- 2p q - (p - q )t)] cos p cosh q 
- [P 4 + q 4 - (P 2 - q 2 ) T|] - 0 

where 

(7.19) 

The second equation in (7.18) is the frequency equation for n » », (no 
Coriolis forces [72]), and the first equation, to the first order of approxi- 
mation in ^/(p Vn) = 0(c), presents the effect of sufficiently small Coriolis 
forces. We note that A^ and ^ are both independent of c and, therefore, we 

may directly seek values of u> and F which make them identically zero. This 
is illustrated in Fig. 7.4 for cr 51 1.5, where the critical load is found to 
be y * 1.67. In Fig. 7.5, the critical load y is plotted against a for suf- 
ficiently small Coriolis forces (the dashed curve). The solid curve for 
torsional flutter in this figure is for the limiting case of n = ® [72]. We 
note that the existence of Coriolis forces does not alter the region of di- 
vergent motion, as is expected. However, it makes this region a closed set, 
that is, in the presence of Coriolis forces, the points on the divergent curve 
indicate neutrally stable states. The horizontal solid line in Fig. 7.5 de- 
notes the limit of transverse flutter for n « », and the horizontal dashed 
line indicates that limit for sufficiently small Coriolis forces [32], (for 

EI x r 2 /C 1 - 1.5). 

It may be of interest to obtain the critical values of y for g * » and 
n « 1. This, of course, provides the upper limit of torsional and transverse 
flutter. The dotted curve in Fig. 7.5 represents this limiting case for 6*1. 
We note that transverse flutter, for g “ ® and n * 1, occurs at y » 47, which 
is not shown in Fig. 7.5. 
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7.3 Stability of a Bar in Parallel Fluid Flow. Taking into Consideration 
the Head Resistance 


The problem of a cantilevered bar placed in a fluid flow was analyzed by 
Kordas [88]. The head resistance was assumed to be represented by pure fol- 
lower force, which was not further related to any parameters of the bar or 
the fluid. Piston theory [89] was assumed to characterize lateral pressure 
on the bar. The continuous system was replaced by a system with two degrees 
of freedom and stability limits in terms of relevant parameters of the problem 
were calculated. 

7.4 Stability of Bars Subjected to Radiant Heat 

In a recent paper Augusti [90] has suggested a special construction of 
the links of an articulated bar, which would make it sensitive to radiant 
heat. The links are made up of cells; the heat absorbed by each cell causes 
thermal deformations which induce a relative rotation of the two adjacent bars. 
An interesting feature of the resulting equations, e.g., for a bar with two 
degrees of freedom, is that derivatives of generalized coordinates (angles of 
rotation) up to third order are introduced. Depending upon the combination 
of relevant parameters, stability can be lost by flutter or by divergence. 

Thermally induced vibration and flutter of flexible booms were discussed 
by Yu [91] and commented on by Augusti [92], where further references to this 
phenomenon can be found. 
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CHAPTER VIII 


LABORATORY EXPERIMENTS AND MODELS 


8.1 Introduction 


As already mentioned, it is the Intent to discuss in this monograph only 
those problems involving follower forces which do not belong to the now almost 
classical areas of aeroelasticity and stability of rotating shafts. If we 
omit these two categories, the only two remaining areas of problems with 
follower forces which have been realized to date involve internal flow through 
flexible pipes and fluid jets impinging on a deformable structure. 

As regards the former area, mention should be made above all of the pio- 
neering work by Benjamin [80,93] on the dynamics of a system of articulated 
cantilevered pipes conveying fluid in which both divergence and flutter were 
observed and stability boundaries were determined analytically and also by 
means of quantitative experiments. Benjamin* s work was continued by Gregory 
and Paidoussis [87,94] who studied theoretically and experimentally continuous 
tubular cantilevers conveying fluid. One should also recall the earlier work 
by Long [95] on vibration of a tube containing flowing fluid, who, however, 
did not observe any instabilities, being interested only in the influence of 
the fluid flow on frequencies of vibration. Divergence of a simply supported 
pipe conveying fluid was observed more recently by Dodds and Runyan [96]. 

Simply supported and cantilevered pipes conveying fluid were investigated also 
by Greenwald and Dugundji [97]. 

The dynamics and stability of slender cylinders surrounded by, rather 
than containing, flowing fluid was studied analytically and experimentally by 
Paidoussis [98-100]. He points out that, provided the flow direction coincides 
with the axis of the cylinder at rest, then, for small motions about the position 
of rest, the forces exerted by the fluid in the two cases of external and in- 
ternal flow are closely similar. This becomes evident on considering that the 
forces exerted by the fluid, excepting those due to fluid friction, in both 
cases arise from lateral acceleration of the flowing fluid, caused by lateral 
motion of the cylinder. In external flow, this acceleration is suffered by 
the virtual or "associated' 1 mass of fluid, which is dynamically equivalent to 
the contained mass of fluid in internal flow. Hawthorne [101], taking advan- 
tage of this similarity, investigated the stability of flexible tubes towed In 
water and demonstrated that divergence instability is possible in such systems. 

In this context it is deemed appropriate to mention some related but con- 
siderably more complex hydro-elastic-pneumatic problems arising in structural 
dynamics of launch vehicles studied by Runyan, Pratt and Pierce [102], as well 
as the broad area of propeller-rotor whirl flutter, a comprehensive review of 
which was recently prepared by Reed [103]. 
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The mechanics of impinging jets was studied to date from the point of 
view of fluid behavior, one of the goals being the determination of the pres- 
sure distribution at a rigid surface. By contrast, in a recent study [104] 
the interest centered on the behavior of a primarily elastic structure sub- 
jected to an impinging jet. Quantitative experiments were carried out and 
compared with theoretical predictions, as described in the following section. 

8.2 Instability of a Mechanical System Induced by an Impinging Fluid Jet 

8.2.1 General 

The mathematical model of the physical system considered here may be 
called Reut's problem, mentioned already in Sect. 6.2. It consists of a can- 
tilever with a rigid plate at its free end, which is normal to the axis. It 
is subjected to a force, acting on the plate, which is always collinear with 
the undeformed axis of the cantilever, Fig. 8.1. Bolotin [8] reports that 
this problem was first posed by Reut in 1939 and solved by Nikolai in the same 
year. In this context, Bolotin suggests that the force in Reut's problem may 
be realized by an impinging jet of absolutely inelastic particles, since the 
kinetic energy of the particles is completely absorbed upon impact. It appears, 
however, that no attempt was ever made to follow up these suggestions, or to 
realize Reut's problem in any other way. Bolotin also suggests that the 
pressure from a jet of liquid or gas may induce such a force when the incli- 
nation of the force, as the bar deforms, is neglected.* 

In an attempt to construct models based on these ideas, it was disco- 
vered that by covering the plate with screens of certain mesh sizes a problem 
very close to the Reut's one may be realized. The resultant force, in this 
case, has an inclination which can be controlled by a suitable arrangement of 
screens of various mesh sizes; the point of application of the resultant force, 
however, always lies on the axis of the undeformed cantilever. When this force 
stays normal to the end plate, the system loses stability by divergence (attain- 
ment of another equilibrium state) ; the force is conservative. On the other 
hand, if the force stays collinear with the undeformed axis of the bar, the 
loss of stability occurs by flutter (oscillations with increasing amplitudes) ; 
the force is nonconservative. By controlling the inclination of the force, 
various degrees of nonconservativeness may be attained. 

The experimental results are obtained using a system with two degrees of 
freedom, rather than a continuous cantilever. The applied force is induced by 
an impinging air jet. The degree of nonconservativeness is controlled by em- 
ploying suitable end attachments, resulting in either divergent or flutter- type 
motions of the system. Also, the effect of viscous damping forces is investi- 
gated. It is found that the experimentally obtained flutter load corresponds 
rather closely to the theoretical prediction when small dissipative forces are 


This, of course, is not acceptable, since it is precisely the presence 
of the component of the force in the direction normal to the impinging fluid 
that, in this case, renders the system conservative. 
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included; this confirms the earlier findings that small damping forces may 
have a destabilizing effect. 

8.2.2 Description of Model and Supporting Equipment 

The model consists of two like rigid rods. Fig. 8.2. One rod is elasti- 
cally hinged to the first rod and free at the other end. The system is con- 
strained to move in a horizontal plane, being supported by long, light wires. 
Various rigid attachments can be placed at the free end of the second rod. 

The attachment consists basically of a rigid flat plate covered with a combi- 
nation of screens of various mesh sizes. This attachment is rigidly fixed 
and mounted normal to the axis of the second rod. In the absence of any dis- 
turbance, the system is in equilibrium when the two rods are collinear (undis- 
turbed configuration) . 

A fixed nozzle is placed along the equilibrium axis of the system, one 
inch away from the attachment, and an air jet is made to impinge upon the 
attachment. The flow rate can be varied by means of a valve. The dynamic 
pressure at the nozzle corresponding to a given flow rate can be read from a 
dial gage. 

It is observed that as the flow rate, and hence the force on the attachment, 
is increased and passes a certain (critical) value, the system does not remain 
in the undisturbed configuration. Stability is lost by either flutter (oscil- 
lations with increasing amplitudes) or by divergence (buckling - the attain- 
ment of another equilibrium state), depending on the nature of the attachment 
used. If the attachment is a flat plate with a smooth surface (a flat sheet 
of aluminum) facing the air jet, then the loss of stability occurs by diver- 
gence. By contrast, flutter- type motion is observed if the attachment is a 
plate with screens of certain mesh sizes placed on the surface that faces the 
impinging fluid. The sequence of photographs in Fig. 8.3 illustrates the 
flutter-type motion, while Fig. 8.4 depicts a buckled state (divergence). Fig. 
8.5 and Table 2 present the numerical values for all the relevant properties 
of the system. 

The supporting equipment consists of a calibrating system which is used 
to correlate the dynamic pressure, hence the flow rate, with the actual force 
which acts on the system. Three square steel plates are placed horizontally 
one above the other, and are separated and supported by sets of steel leaf 
springs. The steel leaf springs connecting the two lower plates permit dis- 
placement in only one direction, while those connecting the upper two plates 
permit displacement only in the perpendicular direction. Two stages are thus 
formed. The displacement of each stage is, with a high degree of accuracy, 
proportional to the component of the force which acts along the direction of 
the displacement. With the aid of a pair of strain gages attached to the steel 
leaf springs, and using a compensating network, readings can be taken which are 
proportional to the respective displacements of each stage. In this manner, 
strain-gage readings can be related to the magnitude of the force acting on 
the system. 
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The supporting equipment just described is used to find the direction and 
the magnitude of the force on the attachment when the dynamic pressure of the 
impinging air jet at the nozzle is known. The attachment is mounted on the 
top plate of the supporting stages and then subjected to the air jet at a 
given angle of incidence, Fig. 8.6. The magnitude and the direction of the 
resultant force corresponding to a given angle of incidence and for a given 
dynamic pressure are thus obtained experimentally. 

8.2.3 Theory 

As already mentioned, the problem of a cantilever with a rigid cross plate 
at its free end and subjected to a force which is always directed along the 
initial, undeformed axis of the cantilever, was first posed by Reut in 1939. 

It is essential to note that the applied force in Reut's problem is not attached 
to a material point of the system, but rather to a line in space. In struc- 
tural mechanics, boundary-value problems are commonly posed for surface trac- 
tions which are connected to the material points upon which they act. As a 
result, the difference between the displacements of the material points and 
of the points of application of the forces disappears. 


In the present problem, the force is induced by the action of an air jet 
upon the end plate. It may be assumed that such an action is equivalent to a 
resultant force whose point of application lies always on the axis of the un- 
deformed system; that is, along the direction of the flow. This force continu- 
ously disengages from the material point on which it is instantaneously acting. 
This force is conservative only if it stays normal to the end plate as the 
system deforms. In the subsequent analysis, we will denote this force by P 
and the angle by which it rotates, as the system deforms, by 

We consider small lateral motions of the system as shown in Fig. 8.5. 

The rigid bar, designated by I, is connected to the support by a rotational 
spring of stiffness and carries at its other end a rotational spring of 
stiffness to which is attached another rigid rod, designated as II. In ad- 
dition, rods I and II are connected to two linear coil springs as shown in 
Fig. 8.5. Since the displacement of the spring connected to bar I is not 
coupled with the motion of bar II, the stiffness properly accounts for the 

effect of this spring. The spring connected to bar II is located at a dis- 
tance d^ from the center of the middle joint and has stiffness K^. 


The inertial properties are represented by seven masses nu , j = 1,2,... 7, 

and seven centroidal moments of inertia I ^ , j = 1,2,. ..7. The mass of the end 

rotational spring is denoted by m^ , and that of the rod I is denoted by m^. 

The central rotational spring has in effect two masses ra^ and which are 

attached to the rods I and II, respectively. The mass of the rod II is m^, 

and m, is that of the collar which fits the attachment having mass m-. 
o / 
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The distance between the centers of the end and the middle rotational 
springs is denoted by j while the mass is at a distance i 2 from the cen- 
ter of the middle joint. The dimensions a^, b^, and are the distances from 
the center of the end joint to masses m^, a^, and a^, respectively, while a^, 
b ^ 9 and c 2 designate the respective distances of m^, m^, and m^ from the cen- 
ter of the middle joint. 

The two rotational springs were made of high tempered spring steel with 
identical geometry and, therefore, they have small damping with, plausibly, 
the same damping constant Since the attachment has a large surface area 

which moves relative to the impinging air jet, an external linear damping with 
constant e 2 appears to be a reasonable representation of the corresponding 

damping mechanism. 


The magnitude of the force due to the impinging air jet is P, the direc- 
tion of which encloses an angle CKp 2 th the undeformed axis, cr is assumed 

to be a constant which will be determined experimentally with the help of the 
auxiliary equipment as described in Sect. 8.2.4. and cp 2 are the respective 

rotations of bars I and II from the initial straight position. 


The following equations of motion are obtained by employing D'Alembert's 
principle: 


^11^1 ^ A 12^2 B 11*^1 + B 12^2 ^ll" P ^ l)^l + (^ 12 * P ^1°^*P2 * ® 

A 21^1 + A 22^2 + B 21^1 + B 22^2 + ^ C 21 + + ^ C 22 + ^ 2*^2 “ 0 
where 

2 2 2 2 

A 11 = ^ m 4 + m 5 + m 6 + ®7^l + m l a l + m 2 b l + m 3 c l + I l + I 2 + I 3 

A 12 " *21 = (m 4 a 2 + m 5 b 2 + m 6 C 2 + ^2^1 

^22 = m 4 a 2 + m 5 b 2 + m 6 c 2 + m 7^2 + I 4 + I 5 + X 6 + I 7 

B 11 * e 2*l 2+ 2 ®1 

B 12 " ®21 = * 2 l l l 2 ' *1 


( 8 . 1 ) 


+ €, 

22 2 2 1 


Hi " K l + V Vl 


( 8 . 2 ) 

cont. 
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( 8 . 2 ) 


C 12 m C 21 * K 2 + K 3^1 d 2 


C 22 " V K 3 d 2 


e l 


Undamped System - Flutter . Consider first the undamped case: i.e. , let 

C 2 = 0. Then =* 0. Assuming solutions of the form 


<P X * $je 


itut 


cp 2 - * 2 e 


loot 


(8.3) 


where i =■ /-l, and are undetermined amplitudes, oj is an undetermined 
frequency and t is the time variable, the associated frequency equation is 

a«) 4 + bw 2 + c » 0 (8.4) 


where 

a = ^1^2“ A 12 


b = 2A i2 C 12 + A 12 W l^ 1+0 ^“ A ll P ^2 a “ A 11 C 22* A 22 C 11 + A 22 P- *1 

c ■ c u n 2 a ' p2 -ei“Cei + ^2 ) ' C 12 W 1“ “ C 12 W 1 + C 11 C 22‘ C 12 2 ' C 22 W 1 


Flutter occurs if u) is complex with a negative imaginary part. The 
threshold (critical) value of P, called P^, is obtained by setting 

b 2 - 4ac - 0 (8.6) 


and is 


P *1 2 = 2 T " £g ± yh 2 k 2 -hkfg-4h 2 jm + hjg 2 + hmf 2 

* f - 4hj f - 4hj y 


(8.7) 


where 


f = A 12 £ 1 (l-h») - A n V + A 22 f L 

8 = 2A 12 C 12 " A 11 C 22 " ^ 2^11 
h - A 11 A 22 - A 12 


( 8 . 8 ) 

cont. 
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j - - l jcr (£ j+ j l 2 ) 


k » C^j^cr - C ^i 01 ~ C 12 jR 1 ” C 22^1 


* “ C 11 C 22 ‘ C 12 


As the value of P is increased, flutter will occur when P becomes equal 
to the lower value of P^. Note that P^ is a function of a through equations 
(8.8)* P* exists only when the argument of the square root in Eq. (8,7) is 
nonnegative. 

Damped System - Flutter . Using an assumed solution of the form (8.3) in 
Eqs. (8.1) results in the following determinant which is set equal to zero for 
a nontrivial solution: 

2 

- U) A^ + + 

- tt)^A 10 + Pi .of + C 10 + iu)B 10 

12 1 12 12 * 0 (8.9) 

2 

- tu A ^2 + c i 2 + + *-^12 

2 

- (jo A 22 + + C 22 + iu)B 22 

If we neglect the product of and $2 * n the expansion of (8.9), we 

obtain two equations by separating the real and imaginary parts. The first 
equation is the same as equation (8.4). The equation resulting from the 
imaginary part yields the following relation: 

2 B 11< M 2^22> + B 22 (C ll" Pi l ) ' B 12^ PX l (1+ “ ) + 2C 12^ ln1 

^22 B 11 + A 11 B 22 ” 2A 12 B 12 

2 

Substituting u) from Eq. (8.10) into Eq. (8.4) and denoting by P^ the thresh- 
old values of P for this case, we obtain 

Y " " h ± kP r ~ 4wv (8 - n) 


where 



(8.12) 

cont. 
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2 +S T* m 


s 


( 8 . 12 ) 


and 

q =* X(3co f - 2) 

r - (1 + 2 c)C 22 + (1 + e)C n - 2C 12 (1 - e) (8.13) 

s = (1 + 2 c)A 2 2 + (1 + e)A 1;L - 2 a 12 (! “ e) 

where 

e s and ^ a » ^2 (8.14) 

Thus, the critical force depends not only on cy, but also on e, essentially 
the ratio of the damping coefficients. The critical force is the lower of the 
two values of P^ and it exists only when the argument of the square root in 

Eq. (8.11) is nonnegative. 

Divergence . For divergence, or buckling, u) is set equal to zero in Eq. 
(8.4). The condition is then 


c = 0 


(8.15) 


Denoting the value of P at which this occurs by P^, we have 


1»2 


= [l 

2 j * 2 j V* 


k^- 4jm 


(8.16) 


where j, k, and m are defined by Eqs. (8.8). 

As are P^ and P^, P^ is also a function of <y, but it is independent of 
the mass distribution. exists only if k^ - 4jm £ 0. 


Results . With the system parameters given, including the spring constants, 
which are determined experimentally (see Sect. 8.2.4), Eqs. (8.7), (8.11), and 
(8.16) must be solved for P for each specified value of a. This repetitious 
task was performed with the aid of a CDC 3400 computer. 

As can be seen in Fig. 8.5, for or = 0, the force P is always directed 
along the equilibrium line; i.e., the line defined by ep^ “ <P 2 = When 0 = 1 , 
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the force is always perpendicular to the surface of the attachment. As dis- 
cussed earlier, in the former case the force is nonconservative, while in the 
latter it is conservative. It turns out that with the present setup, experi- 
mentally realizable or are in the range 0.23 £ a < 0.91. 

Unfortunately, mechanical failure of the joints occurred during the ad- 
vanced stage of experimental measurements and, consequently, when the model 
was reassembled, the spring constants K^, K^* and changed. Thus it became 
necessary to designate the previous model by system'll and the reassembled mo- 
del by system II. With due respect to the difference in system parameters, 
stability curves, P versus a 9 are shown in Figs. 8.12 and 8.13. 

8.2.4 Experimental Procedure and Results 

Correlation of Force With Air Pressure and Determination of <y . To find 
the magnitude and the direction of the force acting on the attachment due to 
a given airflow rate, the supporting equipment described in Sect. 8.2.2 is 
used. 


The nozzle assembly is detached from the model and mounted adjacent to 
the calibrating device. Fig. 8.6, parallel to the direction of motion of one 
of the stages. The rigid attachment is separated from the model and mounted 
on a special bracket on the top plate of the calibrating stages. This bracket 
may be rotated so that the angle between a normal to the attachment and the 
center line of the nozzle, namely 92 * ma y be varied. Markings are provided 
for q>2 = 10 » 15, 20, 25, and 30 deg. 

The first step is to find a relation between the displacement of the 
stages and the force applied to the top plate. This is done by applying known 
forces along the deflections of each stage and noting the strain-gage readings. 
If the direction parallel to the nozzle is designated by x and the perpendi- 
cular direction by y, relations of the form 


P = S.Ae 
x lx 


P = S 0 Ae 
y 2 y 


(8.17) 


may be written. 


P and 
x 


are the differences 


P are the forces, and Ae and Ae 
y x y 

in strain-gage readings between no load and full load, for the x and y-direc- 


tions, respectively. S, and 


S£ are the proportionality constants. 


The next step is to correlate the force, P, with the air pressure, p. 
From the free -body diagram of the attachment mounted on the calibrating sys- 
tem, Fig. 8.7, the following relations are obtained: 


P x ( p ) =* P(p) cos 0f(p)cp 2 
P y (p> = P(p) sin or(p)cp 2 


(8.18) 
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where the force P has been split into its components P and P , which are 

x y 

functions of the pressure, p. The parameter a is assumed to be a function 
of p also. From Eqs. (8.18) we can write 


arc tan 


r P n 

oe «] 


0f(p)«p, 


(8.19) 


For a given attachment and angle strain-gage readings are taken for 

a set of pressures. These in turn yield the forces P x and P^ corresponding 

to each pressure. The angle of incidence, cp^, is then varied from 0-30 deg 

in 5-deg increments and for each value of cp 9 an average value for P /P is 

z. y x 

obtained over a range of pressures p. It turns out experimentally that P^ 
and P^ are linear functions of p, as one would expect, and thus the ratio 
P /P is independent of p. This means that a must be independent of p because 

y * 

of Eq. (8.19). If arctan P /P is plotted versus cp 9 , the result is (very 

y x z 

nearly) a straight line and, therefore, the slope may be interpreted as a in 
Eq. (8.19). a is a constant for a given attachment. 


The critical force is read, or interpolated as the value of at =* 0 
corresponding to the critical value of pressure. For small cp^, P « P x ; this 
is within the scope of the linearized theory. 


In this manner, the value of cr is obtained experimentally for each at- 
tachment. 


Determination of Stiffnesses. 


Dynamic Method. The spring constants K^, » and may be determined 

experimentally by a simple dynamic analysis of various motions of the system. 
The spring constant associated with the end joint and the linear spring 

attached to bar I may be evaluated by locking the middle joint so that the 
two bars move as a rigid unit, Fig. 8.8. After giving a small disturbance, 
the natural frequency is measured from which is determined. In a similar 

manner, the spring constant ^ °f the middle joint may be determined by lock- 
ing the end joint, removing the linear spring attached to bar II and allowing 
the system to oscillate freely. Fig. 8.9. 

Spring constant can be found if is known. The linear spring is 

attached to bar II in its original position and the natural frequency is meas- 
ured. This gives an expression for the combined stiffness from which may 
be evaluated. 
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Static Method, An alternate procedure is to use a static method where- 
by forces are directly applied and the resulting deflection measured. The 
procedure is divided into three steps and is explained in Figs. 8.10 and 8.11. 

Theoretically, these two methods should yield identical results. Experi- 
mentally, the results of the two methods differed slightly. Table 2. The sta- 
tic measurement is to be preferred because the dynamic method depends upon the 
square of experimentally measured frequencies which are not known with great 
accuracy. 

Summary and Results . The basic steps in the experimental procedure are 
as follows: First, choose an attachment and mount it on the model. Raise 

the air pressure slowly from zero and note the critical pressure at which the 
system starts exhibiting amplified oscillations (flutter) or shows a static 
loss of stability (buckling). The supporting equipment is then used to find 
a and to find the force P corresponding to the critical pressure p. The 
spring constants are then determined experimentally for use in the theoretical 
analysis. 

When choosing attachments, it is desirable that they all be of about the 
same weight and that a wide range of a be covered more or less uniformly. A 
wide variety of screens and sandpapers were weighed and combinations were 
chosen that met these requirements. The values of a which were experimentally 
realized lie in the range 0.238 to 0.913, the latter being for an attachment 
consisting of a smooth flat plate. 

In Figs. 8.12 and 8.13, the experimental results are shown together with 
the theoretical curves. As was mentioned, two systems had to be considered 
because of a mechanical failure of the joints. For each experimental run a 
point of instability is drawn on the diagram at the corresponding a and P. 

A 0 is used for a flutter point, while & is used to denote divergence. 

The measurements are labeled 1 through 8 for system I and 1 through 12 for 
system II. 

Table 3 summarizes the experimental and theoretical results and provides 
a comparison between these results. 

8.2.5 Discussion of Results. Conclusions and Recommendations 

The results of this study are summarized in Figs. 8.12 and 8.13 and in 
Table 3. It is noted that the experimentally determined critical points lie 
somewhat below the theoretical stability curves for undamped flutter and di- 
vergence. In the discussion which follows, the possible reasons for this 
discrepancy are explored. 

One of the primary reasons for the discrepancy between the theoretical 
stability curve for undamped flutter and the experimentally observed points 
of flutter appears to lie in the fact that damping is present in the physical 
system. The damping mechanism assumed in the analysis has already been dis- 
cussed. Stability curves for flutter with small damping taken into account 
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are shown in Figs. 8.12 and 8.13 for several values of the damping ratio c. 

It is seen from these figures that in the presence of damping the theoretical 
stability curves come to pass very near the experimental points. No attempt 
was made to determine e with a high degree of accuracy since the assumed 
damping mechanism, while reasonable, was chosen mostly for its simplicity 
and it is doubtful that it represents completely the actual damping in the 
system. Supplementary experiments indicated that the assumed values of e are 
realistic. 

The results presented indicate that damping has a destabilizing effect on 
the system and that the presence of damping extends the flutter region to 
higher values of a. Also, the lower values of the damping ratio are associ- 
ated with lower values of flutter loads and a wider flutter range. This con- 
firms results shown previously in Chapter IV. 

The theoretical curves bounding the regions of flutter (with and without 
damping) and divergence were found to be rather insensitive to small changes 
in system parameters, as indicated in Table 2, with the possible exception of 
the spring constants. The dynamic measurement of the spring constants pro- 
vides another possible source for the discrepancy since the calculation de- 
pends on the square of a measured quantity; i.e., the frequency of free oscil- 
lations. But, the spring constants were determined also using the static 
method previously described. Difficulties may arise here, however, in measur- 
ing the applied force by means of hanging weights on a light string which 
passes over an air bearing. 

Since the two methods of measuring the spring constants gave somewhat 
different results. Table 2, it was decided to investigate the effect of a 5 
percent difference in either K^, K^, or K^. A computer program was written 

in which each calculated spring constant was subjected to a ± 5 percent un- 
certainty. If an envelope is drawn about the nine curves thus obtained, the 
effect is roughly to give a maximum error of ± 6 gm (or ±4-10 percent). 

No other system parameter. Table 2, is subject to an error approaching 5 per- 
cent, except possibly the moments of inertia, but these are insignificant when 
compared to the mass-times-distance- squared terms to which they are added. 

The observed discrepancy between the theoretical curve for divergence and 
the experimental points may be due also, in part, to the uncertainty in the 
values of the spring constants, but the major cause of error seems to lie in 
the possibility of initial imperfections and nonlinear effects. 

Since the physical model is not an ideal linear system free of imper- 
fections, there is no single, sharply defined divergence load. An arbitrary 
criterion of the load required for a one-inch deflection of the middle joint 
was used as the condition for divergence. By this definition, the experimen- 
tal points of divergence were somewhat below (15-25 percent) the divergence 
curves obtained from the linear analysis. Figs. 8.12, 8.13, and Table 3. In 
an attempt to explain this discrepancy it seems advisable to investigate the 
nonlinear divergence theory as well as the effects of initial imperfections. 
This is discussed in detail in the next section for a = 0.717 (run 11). 
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The results of this investigation are shown in Fig. 8.14, with a detailed 
description of the curves given in Sect. 8.2.6. It is noted that the postu- 
lated criterion for divergence gives very nearly the same load for both the 
linear (curve A) and the nonlinear (curve B) cases, and thus the theoretical 
divergence curves given in Figs. 8.12 and 8.13 actually represent the diver- 
gence loads for the nonlinear theory in conjunction with the adopted criterion. 

The strong effect of imperfections on the divergence load is discussed in 
Sect. 8.2.6. Initial imperfections in the amount cp^Q “ 0.01, cp^g ” “ 0.01, as 

shown in curve D, are indeed reasonable for this model. This corresponds to 
a no-load deflection of about 0.1 inch at the middle joint. This small imper- 
fection lowers the theoretical divergence load by about 15 percent. 

Curve F is the experimental force-deflection curve for run 11. Note that 
the shape of the curve differs somewhat from the theoretical curves shown. It 
should be pointed out that the points used to draw this curve are rather dif- 
ficult to obtain since holding the air pressure constant to obtain a deflection 
reading does not prevent the motion of the model. Since the run of the curve 
F is somewhat different from the other curves, the likelihood exists that 
other sources for the discrepancy may be present. It may be appropriate to 
mention here that it has been noted repeatedly in the past that structural 

systems buckle at loads below those theoretically expected. 

To provide better insight into the discrepancy under discussion, the ex- 
perimental procedure was also scrutinized. The method of correlating the air 
pressure as read on the dial gage, to the actual force on the attachment, was 

studied with the conclusion that no appreciable error could be introduced. 

8.2.6 Nonlinear Divergence Analysis 

The equations of motion, assuming cp^ and cp£ are not small, neglecting 

inertial effects, thereby restricting the equations to use for divergence ana- 
lysis, and allowing for imperfections by assuming that the equilibrium confi- 
guration is not a straight line, are 

- k 2 ^2 - 9i> -* p V OS “9 2 Sin 'Pi 

+ K 3 -^ 2.(^1 sin "tPi + d 2 sin ^P2^ cos *9i 


+ sin acp 2 cos = 0 


K 2^2" 9p + P sin Qfc P2 t^2 cos 92 

+ tan sin 9i + ^2 sin 9fc^ + Fj ^l C0S “92 sin 9l 

+ sin 9i + d 2 sin 92^ cos 92 ** 0 


( 8 . 20 ) 
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where - cp x - ep 1Q , q> 2 *= q> 2 - q^ Q , and cp 1Q and cp 2Q are the no-load values 
of cp^ and respectively. 

Restricting the magnitude of cp^ and q^ by setting 


sin cp^ * q>^ - cpi^ /6 

2 

cos c*cp^ * 1 - (cwp 2 ) /2 

the equations may be written as polynomials of the form 

3 3 2 2 2 2 

A l cp l + A 2^ + A 3*l^ + A 4^1 <P2 + A 5*l + A 6^2 

+ A 7 cp lC p 2 + A g cp^ + Agc^ + A 1q « 0 

3 3 2 2 2 2 

B l cp l + *2^2 + B 3 c Pl < ^ + B 4<Pl <P 2 + B 5^1 + *6^ 

+ Bycp^ + Bgq^ + B 9 cp 2 + B 10 - 0 


( 8 . 21 ) 


( 8 . 22 ) 


A computer program was written to solve these two third-degree algebraic 
equations simultaneously for various values of P, aru * ^20 m resu ^ ts 

are given in Fig, 8,14 for or * 0.717 (run 11) in the form P versus q>^. The 
variation of cp 2 with P is essentially similar. The vertical dotted line re- 
presents the angle cp-^ corresponding to one-inch deflection of the middle 
joint, which is the buckling criterion used in this study. 

Curve A represents the linear case for cp^Q 32 q^Q “0. No deflection 

occurs until the buckling load is reached. Curve B represents the imperfection- 
free nonlinear case where the approximations (8.21) are used. The buckling 
loads predicted by curves A and B are rather close. 

Curves C, D, and E are drawn for the values of cp^g and cp 2 g indicated. 

Note that the buckling loads, as determined by the intersection of the re- 
sponse curves with the dotted vertical line, depend significantly on the mag- 
nitude of ep^g and q^g. 

Curve F is the experimental response curve for the model with the attach- 
ment used for run 11 (or * 0.717) in place. 
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8.3 Demonstrational Models 


Considerable insight into the possible types of dynamic behavior of 
mechanical systems subjected to nonconservative forces may be gained not 
only through quantitative experiments, but also by qualitative observations 
of demonstrational models. A set of such models has been recently designed 
and constructed at the Structural Mechanics Laboratory of Northwestern Uni- 
versity [105], and it is intended to develop this set further at the Applied 
Mechanics Laboratory of Stanford University. A brief description of the 
models follows. 

Model A 

The model consists of two like rigid pipe-segments (Fig. 8.15a). The 
first is elastically hinged to a fixed base, while the other is elastically 
hinged to the first and carries a nozzle at the free end. In addition to 
the elastic hinges, the stiffness of the system can be varied by means of 
lateral, spiral springs. The system is constrained to move in a horizontal 
plane, being suspended by long, light strings. A fluid can be conveyed through 
the pipes, entering at the fixed end and leaving through the nozzle. In the 
absence of the fluid, or for small rate of discharge, the pipes are at rest 
and colinear, defining the equilibrium configuration. Two symmetrically pla- 
ced strings in the horizontal plane are attached to the free end of the pipe 
and pulled toward the fixed base at a small angle relative to the pipe axis. 

It is observed that as the flow rate is increased, and passes a certain 
(critical) value, the pipe system does not remain in the undisturbed confi- 
guration. The loss of stability occurs either by divergence or by flutter, 
depending upon the stiffness of the auxiliary coil springs at the free end 
and the tension in the wires. If the coil spring at the free end is suffici- 
ently soft, or is removed, and the tension in the wires small, then the loss 
of stability occurs by flutter-type motion. By contrast, for sufficiently 
stiff coil springs, or for large enough tension in the wires, the system loses 
stability by divergence (Fig. 8.15b). 

In experimenting with this system, it was found that the system can ad- 
mit two distinct critical flutter flow rates. One is associated with rela- 
tively large initial disturbances and the other corresponds to small initial 
perturbations. That is, for a certain range of flow rates, the system is 
asymptotically stable when disturbed by sufficiently small initial input of 
energy, while it oscillates with increasing amplitude about the undeformed 
axis for sufficiently large Initial perturbations (loss of stability in the 
large) . Above this range the system loses stability by flutter for any ini- 
tial disturbances (loss of stability in the small) . 

A thorough and systematic investigation (both analytical and experimental) 
of articulated pipes conveying fluid was presented by Benjamin [80,93]. The 
model described here represents a generalization of Benjamin's system by in- 
cluding a nozzle to control Coriolis forces, lateral springs to control ef- 
fective constraints, and tension wires to control the direction of the resul- 
tant forces acting at the free end. It appears that the existence of loss of 
stability in the large was not observed before in such systems. 
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Model B 


This model consists essentially of a piece of a rubber tube, fixed at 
one end and elastically restrained in the axial direction at the other end, 
at which rotation is prevented. Fig. 8.16. The tube is confined to move in 
the horizontal plane, being suspended by means of long, light strings. A 
fluid can be conveyed through the tube, entering at the fixed end. The other 
end being closed, the fluid is ejected through two nozzles, placed at a cer- 
tain distance from the fixed end symmetrically with respect to the tube in 
the direction parallel to the tangent to the tube at that section. The 
nozzles are mounted in a fixture which is made to slide on an air cushion. 

The sleeve providing the sliding support at the elastically constrained end 
is also supported by an air bearing. In Fig* 8.16 the tubes supplying air 
for the bearings are seen on the left part of the photograph. 

It is observed that the straight equilibrium configuration may be lost 
if the flow rate of the air passing through the tube exceeds a certain cri- 
tical value. Loss of stability can occur by either flutter or divergence, 
depending upon the distance between the nozzles and the fixed end. It may 
be remarked that by attaching a series of pairs of nozzles along the tube, 
the problem of a bar subjected to disturbed tangential follower forces may 
be realized. 

Model C 

This model consists of a cantilevered thin elastic strip at whose free 
end a circular rigid plate is attached in a plane normal to the axis, Fig. 
8.17. The surface of the plate can be varied by placing screens of differ- 
ent mesh sizes. A nozzle whose axis is parallel to the axis of the strip 
can be made to discharge fluid at a constant rate which impinges upon the 
plate. 


It is observed that as a certain critical flow rate is exceeded, the 
cantilever may lose stability by either flutter or divergence, depending 
upon the mesh size of the screen attached to the plate. Both torsional and 
bending deformation are observed to occur for both types of loss of stabi- 
lity, with torsional deformations becoming more pronounced with increased 
eccentricity. 

Model D 


This model consists of a cantilevered thin elastic strip at whose two 
longitudinal edges flexible tubes are attached through one of which fluid 
at constant rate can be conveyed, entering at the fixed end and leaving 
through the open end. Fig. 8.18a. The other tube does not convey any fluid 
and is provided solely to decrease the asymmetry of the cross-section. 

It is observed that as the flow rate exceeds a certain critical value, 
the cantilever loses stability by bending- torsional flutter. Fig. 8.18b. It 
is also observed that a certain range of flow rates restores the original 
undeformed equilibrium configuration which may have been lost by lateral 
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buckling caused by attaching a given weight at the free end. Fig. 8.18c shows 
the buckled configuration at zero flow rate and Fig. 8.18d shows the restored 
original equilibrium position, achieved with a certain flow rate. As the flow 
rate is increased further beyond a certain value, stability is lost by flutter. 

Model E 


This model consists, as in the previous two cases, of a cantilevered 
elastic strip at whose two longitudinal edges flexible tubes are attached. 

A rigid pipe is placed along the transverse free edge and connected to the 
longitudinal tubes. Fig. 8.19. Fluid is conveyed at a constant rate through 
the longitudinal tubes, entering at the fixed end of the cantilever, and is 
discharged through an end opening in the rigid pipe, whose other end is closed. 

It is observed that as the flow rate is increased beyond ascertain cri- 
tical value, stability is lost by bending- torsional flutter. The system may 
be considered as model of an aircraft wing with a jet engine at the free end. 

Model F 


This model consists of a rigid closed cylinder which can roll on a hori- 
zontal plane. A piece of a rigid pipe is attached to the cylinder by means 
of an elastic hinge, which carries a nozzle at the free end, Fig. 8.20a. 

Fluid can be conveyed into the cylinder by means of a flexible tube, which 
then enters the pipe and is discharged through the nozzle. 

It is observed that as the rate of discharge is increased beyond a cer- 
tain value, the system acquires a (stable) equilibrium position such that the 
pipe is vertical and its axis passes through the center of the cylinder. Fig. 
8.20b. As the rate of discharge is increased further, another definite 
(critical) value is reached, beyond which the system begins to execute oscil- 
lations with increasing amplitudes about the preceding equilibrium state 
(flutter) . 

Model G 


This model consists of a rigid cylinder, as in the previous model, which 
can roll on a convex rigid cylindrical segment which in turn is fixed in a 
concave rigid cylindrical segment, this latter being free to roll on a hori- 
zontal plane. Fig. 8.21a. The rigid cylinder is closed at the end planes and 
is provided with an opening and a nozzle on the lateral surface, the axis of 
the nozzle passing through the center of the cylinder. Fluid can be conveyed 
through a flexible tube to the cylinder and is discharged through the nozzle. 

It is observed that as the rate of discharge is increased beyond a certain 
value, the system acquires a (stable) equilibrium configuration such that the 
axis of the nozzle is vertical, Fig. 8.21b. As the rate of discharge is in- 
creased further, another definite (critical) value is reached, beyond which 
the system begins to oscillate with increasing amplitudes about the preceding 
equilibrium state (flutter). Fig. 8.21c. If the convex cylinder segment on 
which the cylinder rolls is replaced by a flat plate. Fig. 8.21d, no flutter 
is observed. 
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Model H 


This model consists of a rigid pipe segment suspended by means of a 
flexible tube and hanging in the vertical position. Fig. 8.22. The lower 
end of the rigid pipe carries an attachment, the essential part of which 
consists of two nozzles placed in a plane normal to the axis of the pipe 
segment, parallel to each other. The flexible tube is connected to a fixed 
base. Fluid can be conveyed through the flexible tube, entering the rigid 
pipe segment and discharging through the nozzles in opposite directions. 

It is observed that for any constant flow rate above a certain minimum 
value, the rigid pipe begins to move like a spherical pendulum with monotoni- 
cally increasing amplitude, which will approach a limiting value for a suf- 
ficiently small flow rate. The minimum value of the constant flow rate which 
produces the onset of the pipe motion is not sharply defined. It is further 
observed that the same motion is initiated if the rigid pipe segment is made 
very short as compared to the flexible tube, and vice versa. 

The problem of a cantilevered bar subjected at the free end to a twist- 
ing moment which rotates with the end cross-section of the bar was first 
considered by Nikolai [2]. He found that the undeformed rectilinear equi- 
librium configuration is unstable for any nonvanishing magnitude of the 
twisting moment. 
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TABLE 2 


SYSTEM DATA 


Dimensions 

(cm) 

Part 

Mass 

(gm) 

Centroidal 
moment of 
inertia 

(fcm-cm^) 

a L « 0.692 

1 

10.20 

rJ 0 

b x * 16.3 

2 

22.0 

1655 

c 1 = 31.9 

3 

42.1 

0 

a 2 = 0.692 

4 

10.20 

- 0 

b 2 = 16.3 

5 

22.0 

1655 

c 2 = 32.3 

6 

3.2 

r+J 0 

d 2 = 25.3 

7 

43.5 

111 


JL ± = 32.4 
i 2 = 32.6 


Spring constants 


System I 


System II 


K. 


Dynamic method 

5.70 x 10^ dyne -cm 
9.12 x 10 5 dyne* cm 
3.50 x 10^ dyne /cm 

5.34 x 10^ dyne* cm 
9.02 x 10^ dyne -cm 

3.35 x 10^ dyne /cm 


Static method 


... 


5.45 X 10^ gm*cm 
9.41 x 10^ gm-cra 
3.53 x 10^ gm/cm 



TABLE 3 


SUMMARY OF NUMERICAL RESULTS 


SYSTEM I SYSTEM II 




Experi- 

mental 

Theoretical 

X 

error 

Experi- 

mental 

Theoretical 

X 

error 

Run 

# 

a 

* 

P 1 2 3 

crit 

(g°0 (gn) (ga) (gn) 


* 

P 1 2 3 

crit 

(gm) (gm) (gm) (gm) 


1 

0.343 

F 

56.4 

72 

55 

• • • 

+ 2.5 

F 

70.2 

89 

62 

... 

+12.9 

2 

0.327 

F 

55.2 

70 

55 

# # a 

+ 0.3 

F 

69.7 

88 

62 

... 

+11.1 

3 

0.560 

B 

94.9 

• * 

* • 

124 

-23.4 

B 

118.3 

. . 

» • 

140 
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Observed loss of stability: F = flutter, B = buckling. 

1 Undamped flutter. 

2 Damped flutter, e *= 5.0. 

3 Buckling. 
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Fig. 1.2 Equilibrium curves of a centrally loaded column 
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Fig. 3.2 Load versus 
of parametei 




Fig. 3.3 Detail of load versus frequency curve for a m 0.5, 
illustrating multiple ranges of stability and 
instability 


159 



Fig. 3.4 Critical loads versus parameter or 
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Fig. 3.5 Regions of stability for > 0 
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Fig. 3.10 Possibilities of flutter end divergence 
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Fig. 3.12 Details of loss of stability 

in areas I through XII (Cont'd.) 
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Fig* 3*12 Details of loss of stability 

in areas 1 through XII (Concluded) 
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Fig. 3.14 Loss of stability in regions I through IV 
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Fig. 4.2 Orthographic projection* and perspective of root curves 
of characteristic equation with no daaping 






Fig. 4.4 Critical load veratia ratio of damping coafficienta for « 1 
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Fig. 4.9 General nature of roots of characteristic equation 
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Fig. 4 .16 Appropriate values of 0 versus values of or for complete elimination 
of destabilizing effect 
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Fig. 4.24 System with two degrees of freedom 
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Fig. 4.25 System with distributed external damping 
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Fig. A 28 Critical force as a function of a for various values of B.. 
Thick line joins values of transition force 
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Fig, 6.1 Stability region in the parameter plane for Example 2 










Fig* 7.2 Types of instability as a function of system geometry 
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Fig. 8.4 Buckled state: 


Divergence 
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Fig. 8.9 Configuration to find and K, by dynamic 
method 
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Fig. 8.12 Stability diagram - System I 
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Fig. 8.16 Demonstration Model B 
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Fig. 8*18 Demonstration Model D 
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Fig. 8.18 Demonstration Model D 
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Fig. 8.22 Demonstration Model H 
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